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I. INTRO D0CTI38 



A. GENERAL PRESENTATION OF THE PHENOMENA 

Great progress has been made in the static and dynamic 
analysis of complex structures through the continued devel- 
opment of discrete element methods of structural analysis. 
Tremendous improvements in computing power have made 
possible the study of fully nonlinear problems. 

The analysis of the response of a structure submerged in 
a fluid, is severely complicated by the intrusion of signif- 
icant f luid-struct ure interaction affects. Recently, the 
development of a variety of surface interaction approxima- 
tions has provided the means for a more efficient analysis 
of the coupling between the structure and the surrounding 
fluid. 

Computer codes for structural analysis are well- 
developed so that the fluid-structure interaction is, for 
the most part, handled by adding new capabilities to 
existing strucrural analysis programs. 

It is a well known fact that the primary threat to a 
submerged structure is the shock wave that results from an 
underwater explosion. However, the complexity of the phys- 
ical phenomena involved, along with the difficulties encoun- 
tered in obtaining experimental results have been serious 
drawbacks for the analysis of these types of problems. But 
there is a definite need for investigations of large defor- 
mations and buckling problems in a structure submitted to an 
underwater explosion. 



B. OBJECTIVES 



This study deals with the nonlinear response of a 
submerged cylindrical shell to a shock wave. The existing 
finite element code EPSA (Elasto- Plastic Shell Analysis) 
[Ref. 1] which includes nonlinear affects and a surface 
interaction approximation was selected for the study. In 
order to alleviate the tedious interpretation of results at 
points throughout the shell, PATRAN-G, a color graphics 
system, was used. P ATRAN-3 allows for a global representa- 
tion of a quantity distribution over the structure rather 
than the discrete representation given by a computer output. 
The objectives of this study were to merge the finite 
element code EPSA with the color graphics system PATRAN-G, 
and to conduct an analysis of the response of various types 
of cylindrical shells to a spherical shock wave generated by 
an underwater explosion. 
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II. THE EPSA COMPOTES PROGRAM 

A. PRESENT ATIOJf 

EPSA (Elas t o-Plastic Shell Analysis) [Ref. 1] is a 
computer program developed by Haidlinger Associates and 
funded by DNA/N AVSEA/ONR for the purpose of the analysis of 
submerged stiffened shells under dynamic loadings. It incor- 
porates a number of specific features which are geared for a 
more efficient analysis. 

In particular, EPSA allows: 

-The analysis of shells in an acoustic medium, subjected to 
both low and high frequency shock loadings. 

-An efficient modelling of the elasto-plastic behavior 
-The inclusion of large displacement effects to analyze 
dynamic buckling situations and post- buckling behavior. 

-The modelling of stiffeners and internal structures. 

-The fluid-structure interaction effect 

The following sections describe the equations of motion 
for a submerged structural shell in an infinite fluid 
sujected to a pressure loading and investigate the modelling 
of the surrounding fluid. The last sections are devoted to 
the finite element discretization as well as to specific 
features concerning EPSA. 

B. EQUATIONS OF MOTION FOR THE SHELL 

Considering a thin shell of thickness h , volume V, area 
R submerged in an infinite fluid (figure 2. 1) . The shell 
stress resultants are defined from the stress tensor by : 
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Figure 2.1 Shell Stress State, 
principle of virtual work gives: 




(s} T { fie} dH 




{p} T (fiu}dR - 




P (uHfiujdR 



( 2 . 1 ) 



where 

{u}=(u 1 ,u 2 ,w ) T is the displacement vector at each point 

{s} = (Nh,N 2 2 » N 12 , M 1 1 ,M 2 2 /Mi 2 ) T is the stress resultant 

vector 

{el®( e n»e 2 2 ,k 22 /2k i2 )T is the strain resultant 

vector 

p is the mass per unit area for the shall. 

The symbol fi will refer to a virtual quantity, and the dot 
or star denotes a differentiation with respect to time. 
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The first tern of equation (2. 1) represents the virtual work 
of internal forces , the second represents the virtual work 
of external forces (i.e. pressure, point loading, etc) , and 
the third one represents the contribution of inertia forces 
in the virtual work. Thus, this last term expresses the 
effects of dynamic phenomena on the structure. 

C. FLUID HODELING 

In the case of a submerged structure, the external 
forces are the pressures applied at the fluid-structure 
interface. 

as the shock wave hits the stucture it gets reflected, 
thus inducing a pressure term p r . In addition, the motions 
of the shell will also generate a radiated wave, with a 
pressure contribution p rS 

Therefore, the pressure at the fluid structure interface is 
the sum of the incident, reflected and radiated pressures: 



P = Pi + P r + Pra 



( 2 . 2 ) 



Where 

p = Total dynamic pressure. 

p i = Pressure associated with incident free field pressure 
wave . 

p r = Reflected pressure due to the interaction of the inci- 
dent wave with the structure, the structure being fixed and 
rigid. 

Pra = Radiated pressure due to the normal movements of the 

surface of the structure in the fluid 

p s = p r + p ra is called the scattered pressure. 

The methods for getting the scattered pressure p s will 
now be investigated. Assuming a spherical wave in an 
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acoustic medium with a sound speed c , the pressure is 
determined by the well known wave equation: 



and the proper conditions at the boundaries of the fluid 
domain . 

One alternative for solving the previous problem would 
be simply to use a finite element discretization of part of 
the fluid domain , imposing a radiation condition at the 
boundary [Ref. 2}. However, this would require a very 
significant fraction of the computing effort that could not 
be devoted to the structural modelling. 

Therefore, a more efficient way of computing the scattered 
pressure is required. 

The Doubly Asymptotic Approximation (DAA) imparts upon 
the structural model a surface loading composed of incident 
and scattered waves. 

In the high frequency limit, it can be shown that the 
scattered nodal force vector (F s } is related to the wave 
particle velocities normal to the structure's surface by 

[Ref. 3] : 



Where [ A ] is the diagonal matrix of nodal areas (areas asso- 
ciated with each node) and [U s } is the vector of nodal scat- 
tered normal velocities. Therefors, in this high frequency 
case the shock wave is simply a plane wave and equation 
(2.4) states that the pressure is proportional to the fluid 
velocity. 




(2.3) 



{F s } = [A]{U S ] 



(2.4) 
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In tha low frequency limit tha fluid structure interac- 
tion is governed by the relation: 

{F s } = [Jl v ] {U*} (2.5) 

£ 

Where {U s }=— {0 S } - s no ^ a ^ normal acceleration vector 

and [ ] is the added mass matrix computed in an incompres- 
sible fluid. 

Thus, in the low frequency case, the loading is due to 
tha motion of the rigid structure in an incompressible 
fluid, a problem typically found in hydrodynamics. 

When a broader range of frequencies is considered, one 
can combine the two previous equations with the differential 
equation governing the scattered pressure [Ref. 3], giving: 



[A}» {F s *} * pcC^Ji {F s } = pCfUs*) (2.6) 

Where {F s *} = A (F s } 

Defining the vector cf nodal scattered pressures (p s J by: 

{p s } = U]* 1 (P s } 

we get: 

[M V ]{P S } + pC[A]{p s } = pc[M v J{Us} (2.7) 

which is the set of equations that gives the scattered pres- 
sure at each node of the wetted surface of the shell. 
Equation (2.7) gives exact results in both the high and low 
frequency limits. Thus, DAA allows the modelling of the 
fluid-structure interaction via a coupled set of differen- 
tial equations at the wetted nodes of the structure instead 
of modelling the whole fluid with a finite element grid. 
The use of DAA will free some memory space in the computer 
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for a better modelling of structural behavior while at the 
same time giving some reasonably accurate results for the 
loading of the structure. 

Therefore, the equations for the study of underwater 
shocks will consist of a coupled set of structural equations 
that come from the Principle of Virtual Work and of fluid 
equations at the wetted nodes that :oie from DA A equations. 



D. FINITE ELEMENT PROCEDURE 



1 • Discre t iza t ion 

The principle of virtual work is rewritten for the 
structure introduced in section B [Ref. 4 ] 



/, 



{srCSeJdR-/ {p} 
R •'R 



- {pf {Su}dR +/ p{uf{6u}dR = 0 

Jr Jr 



( 2 . 8 ) 



The surface of the region is covered by a quadrilat- 
eral mesh of N elements, each having an area A.^ . The 
previous integral can then be replaced by a summation of 
integrals over A^ . : 



f t t 

/ (s) {<5e}dR = ) {s) i {<$e}. dR 

-R i = l A • 



(2.9) 



The integrations over A . ar e then performed by dividing the 

1 

element into four regions (figure 2.2) We have then: 



f Mi { 5 e}. dR = Y, t s k f 5e k>i A 

J a . k= 1 



( 2 . 10 ) 



and therefore equation (2.9) becomes : 

f isf (6e}dS = ^ ^ (s)i {5e k }. A* 

J R i= 1 k= 1 

Using the same procedure, it can also be derived : 



( 2 . 11 ) 
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( 2 . 12 ) 



/, 



CP} T {Su} 



dR= zL Z {P)T C<SUl i A i = (P)T (6q) 

i=l k=l 



l 



p{u) (<Su) d R= 



N 4 




i = 1 k=l 



(2. 13) 



Where [ M ] is the mass matrix , {P} is the vector cf exter- 

nally applied forces, and {q} is the vector of nodal normal 
displacements for the structure. 




Figure 2.2 Grid Points in EPSA. 

By definition, finite element discretization can 
express the displacement {u} at any point of an element as a 
function of the displacements at the corner points of the 
element, defined by {q} . 

{U} = [H ] {q} (2. 14) 

Where [H] is a 6x*2 matrix of interpolation functions. 

Combining the derivatives of {u} will give the strain vector 
{e} . In matrix form, equation (2.14) gives : 
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[B] (q) 



(2. 15) 



( 9 ) = 



Where [B] is a 6^12 matrix function of the geometry of the 
element as well as a function of the previous deformations 
in the case of large displacements. 

Using the previous result, equation (2.9) is rewritten in 
the following way : 



f {s} T {<S e} dR= £ £ (s)^ [B k ]i {6q}i {F} T {6q} 

i=lk=l 

Where {F} is the internal force vector. 



(2.16) 



Combining the previous equations in equation (2.8) 
principle of virtual work becomes: 



CM] fq) 




{FJi I 



, the 



(2. 17) 



Therefore, the principle of virtual work has been tranformed 
into a system of ordinary differential equations which are 
much more amenable to numerical treatment. 

In EPSA, each arbitrarily shaped quadrilateral 
element is defined by four corner nodes, each having three 
translational and no rotational degree of freedom. In order 
to represent the behavior in bending eight nodes not contig- 
uous with the element are also used (figure 2. 3) . Every 
element accesses twelve nodes and has twenty degrees of 
freedom: three translational degrees at the corner nodes 

and one corresponding to the displacement normal to rhe 
surface for each of the eight exterior nodes. The bending 
behavior (second derivative term) is expressed in terms of 
the nodal displacements via a finite difference technique 
[Ref. 4]. 
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Then the nodal displacement vector of element i is 

simpl y: 



(q} = < u 1 ' d 2 ' U 3 ' a 4 ' V l'*‘ ,v 5 ,...v 12 f <2.18) 

Conventional finite element codes utilize three 
translational and two rotational degrees of freedom at each 
node (each element has 5 4 = 20 d.o.f. as in EPS A ) . However 
the masses associated with rotational degrees of freedom are 
very small, leading to numerical instability. The use of 
the aforementioned formulation alleviates this problem since 
only translations are considered and, in addition, the 
number of unknowns is reduced, leading to simpler and more 
efficient computations. 

It must be pointed out that in order to model the 
bending behavior at the edges of the shell, a set of artifi- 
cial nodes has to be created. The finite element grid will 
then consist of the nodes defined in the input file plus the 
artificial nodes along the boundary of the sheet 
(figure 2.3). 

2. Strain Disp la cemen t Relation 

The principle of virtual work deals only with 
strains. Since the finite element approximation gives the 
displacement at each point, the equations relating the 
strains to the displacements are needed. In this study, the 
Donn ell- Vlasov nonlinear kinematic equations of shell theory 
are employed, and the strain-displacement relations are 
described in table I. 

Using equation (2.15) in the equations detailed in 
the previous section will give the finite element approxima- 
tion of strains in terms of nodal displacements. 
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TABLE I 

Donne 11- Vlasov Shell Equations 
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3. Shell C onstitut ive Rel atio ns 

The shell constitutive relations relate the stress 
resultant rate vector to the shell strain rate vector. In 
matrix terms: 
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(s) = [D](s} 



(2. 19) 



Where [D] is the Elasro-Plastic tangent stiffness matrix. 




Figure 2.3 Nodal Points Organization. 

The stress- strain relation used in EPSA differs from 
the classical Elasto-Plastic theory in that the formulation 
involves shell stress resultants rather than stresses at 
points throughout the thickness of the shell. In other 
words, EPSA uses relation 2.19 integrated over the thickness 
of the shell. Thus, there is no need to compute the stresses 
throughout the shell, which results in a significant savings 
in storage space and processing time. However, the stress 
resultants N^and cf the shell theory are not sufficient 
to describe the state of stress, and certain higher order 
moments must be combined with the stress resultants to form 
the dynamic variables of the problem. The coefficients for 
the integrated constitutive equation have been determined 
using experimental results [Ref. 5]. 
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These constitutive equations consist of a yield condition, a 
strain hardening law and a flow rule: 

-The yield condition indicates whether part of the shell has 
started to yield. (figure 2.4) 

-The strain hardening law gives the evolution of stresses in 
the shell after plasticity is reached. 

-The flow rule gives the plastic strain rate in the shell 
after plasticity. 
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Figure 2.4 Yield Situation in the Shell. 



4 . Soluti on Pr oced ure 

EPSA uses an explicit central difference scheme to 
solve the virtual work and fluid loading equations detailed 
in section B. As indicated in appendix B, the explicit time 
integration procedure requires a small time step and is only 
conditionally stable. However when stable, it always 

converges to the exact solution, as opposed to implicit 
schemes that are unconditionally stable but may lead to 
unrealistic results. Furthermore, in problems involving the 
treatment of shocks, accuracy requirements preclude the use 
of large time steps. The central difference scheme seems, 
therefore, particularly optimal. 
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Assuming we know the solution at tine t , the central differ- 
ence scheme applied to equation (2.17) gives the solution at 
time t+A t : 



t-tAt t 2. 

{V}. = {7} + At ( CPI i - £ CF>k ) (2.20) 

1 1 TT k=1 
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[Hef. 6]. 

E. EPSA CAPABILITIES 

The structure to be analyzed is conceptually divided 
into constitutive parts called ’’sheets." Each sheet is a 
curved section of a shell with an arbitrary number of nodes 
and elements (figure 2.5) Its shape is limited to a surface 
that can be described by a smooth continuous function 
without discontinuities in its slope. The elements within 
the sheet are limited to a rectangular organisation (figure 
2.5) . 

Thus, a cylinder with end caps would consist of three 
sheets: a cylindrical sheet and a sheet for each end 

(figure 2.5). Three sheets are required because of the 
slope discontinuity at the edge between the cylinder and the 
end caps. 

Dividing the structure into sheets isolates the diffi- 
culties associated with the boundaries into a set of artifi- 
cial nodes along the perimeter of the sheet. When several 
sheets are required ,EPSA takes care of the compatibilities 
of displacements, rotations and moments along the edges. 
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Figure 2.5 Sheet Organization. 

Thus, any arbitrarily- sha ped structure can be analyzed 
using EPSA by dividing it into a number of sheets. 

Two options for choosing elements are available in EPSA. 
The first option exists to employ a generalized quadrilat- 
eral element. The second option exists to employ a rectan- 
gular element and uses a specialized formulation for this 
type of element. 

The coordinates which can be either cartesian or cylin- 
drical always lie within the plane of the sheet. The z 
direction lies in a positively outward direction normal to 
the sheet. Each sheet contains its own local coordinate 
system, there is no global coordinate system when multi 
sheets are merged (figure 2.6). 

Prior to generating a finite element mesh for a sheet 
one must establish the side numbers of the sheet. The side 
numbering scheme is arbitrary as to the choice of sheet 
number one. However the specification of sides 1 to 4 must 
proceed in a counterclockwise direction when the sheet is 
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Figure 2.6 Coordinate System. 

viewed from the positive z direction. At the intersection 
of side 1 and side 4, element 1 is first defined. Then the 
node/element number is incremented by one until it reaches 
side 2. Then it returns to side 4 and continues the count 
for the next row of nodes/elements. 

Thanks to the exclusive use of quadrilateral elements 
and to the specific numbering procedure, the table of 
connectivity is implicitly defined when generating the nodal 
points, thus simplifying considerably the input require- 
ment. 

The inclusion of structures internal to the cylindrical 
sheet is enacted in EPSA through the use of internal sheets. 
Structures internal to a cylinder are therefore modelled as 
individual sheets having the same capabilities as any 
general EPSA sheet. 

Two types of internal structures are available: 

-Sheets (beams, plates or shells) whose connection to the 
cylindrical shell lies parallel to the axis of the cylinder. 
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-Sheets (beams, plates, shells) whose connection lies in the 
circumferential direction of the cylinder. 

Furthermore, in order to model internal equipment 
(machinery, etc ) which does not deform but contributes to 
the inertia of the system, concentrated masses can be input 
at nodal points. 

The user must be aware that the previously discussed DAA 
is only implemented for cylindrical stuctures. Prior to the 
finite element analysis the user must compute the added mass 
(virtual mass) matrix defined in equation (2.4). This is 
done by using the ACCESION program, which creates a virtual 
mass file that EPSA reads when computing the flu id- structure 
interaction. 

Finally, EPSA in its latest version takes l-ocal cavita- 
tion into account. When the total pressure computed by EPSA 
is found to be negative, it is simply set to zero since 
water cannot withstand any tension. 

F. USING EPSA 

The input file for EPSA can be generated either 
directly, or via the interactive program PEEPSA that prompts 
the user to give the input data . When the input file is 
created, all the data are in free format. 

The nodal points can be generated semi-automatically 
(see the user's manual), and this option is very helpful and 
time saving when generating big models. EPSA can be run 

interactively for small models or on batch for bigger jobs. 
For instance, a cylindrical shell with 1440 elements and 
1517 nodal points takes 0.0129 sec. of VAX CPU per time step 
per element. The whole model would take about half an hour 
for 100 steps. 

EPSA creates an output file in which all input data is 
echoed, and outputs the pressures, stresses, strains, veloc- 



28 



ities, displacements at the nodal points requested by the 
user in the input deck. 

The value of the time increment At should be selected so 
that ; 

At < 1/2 S m . n ( P/E) 55 (2.21) 

Where<5 .is the smallest disxance between the midpoints of 
opposite sides of an element, for all elements of the sheet 
(figure 2.5), and 1/2 is a safety factor. In other words, 
the time step increment has to be less than the time taken 
by a wave to propagate from an element to another. In equa- 

h 

tion 2.21, (E/p) is simply the wave speed in the material. 

The virtual mass array ( VM A) is created on unit 20, 

therefore one should not use this unit for any other purpose 
than HEAD operations. 

Because of the simple organization of its input file, 
EPSA has been found fairly easy to use. The user can perform 
major changes in the model very quickly and efficiently. The 
ACES I ON module has been found to work well except for cylin- 
ders of large dimensions (L=1400' 1 , D=240") for which the 

size of the virtual mass array grows une xpected ely from a 
reasonable 4 blocks to 190 blocks of VAX memory size. 

However, EPSA has been designed for some specific types 
of f luid-stucture interaction analysis and its limitations 
should be pointed out: 

- Only beam type stiffeners can be considered 

The fluid structure interaction is only enacted for a 
circular cylindrical sheet, which takes away much of the 
flexibility the program. 



29 



III. BPSA/ P ATR AN;G INTERFACE 



A. INTRODOCTIOH TO COLOR GRAPHICS SYSTEHS 

In dealing with the response of a structure to a 
loading, quantities such as stresses, strains, velocities 
and displacements must be analyzed at a number of nodal 
points, which makes the interpretation of computer outputs 
very tedious. Color graphics systems offer an effective 
solution to this problem by providing a global representa- 
tion of a quantity distribution over a stucture rather than 
a discrete representation given by a computer output. A 
color graphics system consists of an interface between the 
computer, the user and the color terminal. A typical system 
would be a software package that allows the user to create 
models on the screen as wall as to display any data in terms 
of color intensity. It is known that a picture is worth 
several hundreds of words. Therefore, merging a finite 
element program with a color graphics system would be a 
major improvement in engineering analysis. 

PATRAN-G [Ref. 7] is a color graphics system specifi- 
cally designed for finite elements, it permits the engineer 
and the computer to work together towards the creation of a 
model. The designer creates an image on the screen and 
PATRAN automatically translates the physical model into a 
standard finite element input deck. Another important 
feature of PATRAN-G is its ability to assist the user in the 
interpretation of results through its post-processing facil- 
ities which include: the deformed geometry with magnified 

deformations, the color coding of elements based upon any 
response parameters such as displacements, stresses and 
strains. In particular, the contour levels of the 
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aforementioned quantities can be superimposed on a 

3-dimensional image of the model, thus allowing for a global 
analysis of a complex structure. 

B. HERGIIG EPS A AND PATRAN-G 

As described in chapter I , the structures under study 
have fairly short EPSA input files. Furthermore, in this 
study dealing with fluid-structure interactions on a 

circular cylinder, only structures that consist of one sheet 
are considered. Therefore, because of the simplicity of 
both the input file and the stucture under study, there was 
no need to design a module converting a PATRAN-G model that 
is created on the screen into an EPSA input file. 

The remaining tasks were the following: 

-Display the original finite element model defined in the 
EPSA input file on the screen (original geometry) . 

-Display the nodal points and element results that are 
computed from EPSA on the screen (postprocessing) . 

1 . Origin a 1 G eomet ry 

The input of a finite element model into PATRAN-G 
can be done by creating a "neutral" file, as described in 

PATRAN-G terminology. The neutral file 1 is intended to 
provide a simple link between PATRAN-G and the outside 
world. It is written entirely in 80 character card images 
and all the data is organized in small "packets" of two or 
more card images. Each packet contains the data for a 
fundamental unit of the model such as node coordinates or 
elements definition. Since our only purpose was to display 
the original geometry of the structure, a limited number of 



1 Additional information about neutral files can be found 
in the PATRAN-G user's manual [Ref. 7] 
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data packers (4) w as to be created: 

-File title (packet 25) : two cards containing the user 

title. 

-Summary data (packet 26) : two cards containing the number 

of nodes, elements and the date and time at which the 
neutral file was created. 

-Node data (packet 1): contains all information concerning 

nodes: node number and coordinates in a global coordinate 

frame. 

-Element data (packet 2) : contains the connectivity table 

for the finite element model. 

-End of file (packet 99) : end-of-file cards. 

He have seen in the presentation of EPSA in chapter I that 
the nodes are defined in a local, sheet-attached coordinate 
system, that artificial nodes are created along the edges of 
the sheet to model the bending behavior, and that no connec- 
tivity table was input. Therefore, the translator module 
that would be created had to: 

-skip the set of artificial nodes and properly renumber the 
grid 

-perform a change of coordinates for all local data 
-generate the connectivity table. 

It was decided to employ a modular design in which 
each routine would perform a specific task. A modular design 
would allow further changes to be made quickly and effi- 
ciently by modifying only the relevant routines. The imple- 
mentation in EPSA was made using a series of "calls”, thus 
minimizing the risk of interference with the finite element 
computations. 

The translator module created was made of four 
routines : 
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-PRELIH : computes the number of nodal points, the number 
of elements and displays the first two data packets (25,26) 

-SHEETF : scans through the nodes, skips artificial 

nodes, renumbers the grid, performs the required changes of 
coordinates and displays the node data packet (31) 

-SHCONN : scans through the nodes, connecting each node 

to the elements it belongs to and displays the element data 
packet (02) on the neutral file. 

-ENDNEU : displays the end-of-file packet (99) 



2 • Osina the Translator Module 

The translator calls were implemented in the routine 
REPORT of EPSA . Any run of EPSA creates a neutral file on 
unit 19. The neutral file name is therefore FOR019.DAT if no 
"ASSIGN" statement has been issued prior to the computer 
run. The finite element model might then be displayed on the 
graphic terminal (Ramtek, Tektronix) via the neutral input 
mode of EPSA (see [Ref. 7] for more details). An example of 
finite element model output on Tektronix 4014 is given on 
figure 3.1 . 



>lem e n ta ticn of P o st pr ocessing C apabi lities 



Postprocessing capabilities exist to assist the user 
in the interpretation of computer results. It is simply a 
process of generating displays and reports based upon a 
combination of geometry and the results of an analysis. 

The results cf analyses are transmitted to PATRAN-G 
for postprocessing via "neutral results files" as decribed 
in PATRAN-G terminology. Unlike the model neutral file 
described in the previous section, results files are in 
binary rather than in card image form. 
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Figure 3.1 



Example of 



Finite Element Hodel Display. 



One can distinguish between two major kinds of post- 
processing displays: deformed geometries and element 
guan tities. 

a. Deformed Geometry 



PATRAN-G 
skip the 
we 11 as 
results 
created 



A displacement 
had to be created, 
artificial nodes, 
to write the nod 
file. The disp 
in module NEUDISP. 



results data file required by 
Again, the module created had to 
perform changes of coordinates as 
al diplacements in the neutral 
lacement results data file was 
Its organization is given on 



table II. 



A small module PLOTDISP that decides at which 
time steps the results should be output was created. The 
"call" to PLDTDISP was implemented in module COMPOTE of 
EPSA. The overall structure of the translator module is 
presented on table IV at the end of the chapter. 
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TABLE II 

Organisation of Displacement Results File 






C olumn 



Conten t 



1 X displacement in global coordinates 

2 Y displacement in global coordinates 

3 z " " " 

4 Displacement normal to the shell without 
rigid body mode. 

5 Velocity normal to the shell 



i 



j 



A sample of deformed geometry output on Tektronix 4014 is 
given in figure 3.2. 

b. Element Quantities Postprocessing 

Any element related quantity can be the target 
of postprocessing. In general these types of quantities are 
the results of finite element computations supplied to 
PATRAN-G through the neutral element results file. The 
neutral element results file is different from the neutral 
displacements results file detailed in the previous section, 
however it shares a similar format 'Ref. 7], Element quan- 
tities such as stresses in local and global coordinates and 
von fiises stresses are computed in module NEUSTRE whose 
overall structure is similar to NEUDISP described earlier. 
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Figure 3.2 Deformed Geometry Output. 



The organisation of . the neutral element results file is 
given on table III. 

As described in chapter I, EPSA does not compute 
the stresses through the thickness of the shell. Instead one 
uses the integrated quantities of the shell theory [Bef. 8] 



^ V2 

N i j= I & i j dt and 

Jt 1/2 



h/2 

M i j = / a i j 

Jr/2 



tdt 



One can expect the stresses on the shell to be 
maximum at the extreme fibers. NEUSIRE computes the von 
Mises stresses at the top and bottom fibers and writes the 
maximum value in the neutral file. At the surface of the 
shell no shear stresses are involved. Assuming a linear 
distribution of bending stresses and a uniform distibution 
of membrane stresses, one can easily derive : 
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»ij /h * Mij 6/h 2 



(3.1) 



'3-D 



The first term of the previous equation is the membrane 
force contribution, the second is the bending moment contri- 
bution. The von Mises stresses are then computed using the 
well-known relation: 

°vm = V ( °1 - a l°2 *°2 ) (3-2) 

In a similar way to the displacements results 
file, a module PLOTSTRE that decides whether or not to 
output the element results was created and called from 
COMPUTE. The overall structure of the translator is given on 
table IV. 



TABLE III 

Organisation of the Neutral Element Results File 



Column 


Conte nt 


Description 






22 


st re, 1 


Element local 


stress. 


direction i 


23 


stre , 2 


i« i« 


ii 


" 3 


25 


stre, 4 


Element global 


stress 


, direction x 


26 


stre, 5 


II l« 


H 


" y 


27 


stre, 6 


Element global 


shear. 


direction xy 


31 


von 


von Mises stress 
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c. Displaying the Results Quantities on the Screen. 

The EPS A input deck has been modified so as to 
create the results files (displacema nts, elements) required 
by PATRAN-G. At the end of the second input card the user 
specifies the number of displacement results and the number 
of element results files to be created (at least one) . 
Obviously, those two inputs are also in free format. The 
neutral results files (elements, displacements) will be 
generated at equal time intervals as the computation 
proceeds. The results corresponding to the last time step 
are always output. 

The element results file is created on unit 16 
and the displacement results file on unit 18, thus corre- 
sponding to files FOR016.DAT and FDR018.DAT respectively. A 
new version of those files is created each time an output is 
requested. 

Example : 

If five (5) neutral element results files are requested on 
the input card of a run of 20 steps, five files FOR016.DAT;1 
to FOR0 16. DAT; 5 will be created, corresponding to time steps 
4 to 20 respectively. 

For the displaying of element and nodal points 
results, the user will refer to [Ref. 7] section 20. The 
title of the run (first card of EPSA input deck) will be 
displayed on the screen along with the time at which the 
results have been requested . 
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TABLE IV 

Structure of the Translator 



V 



CALL PL3TDISP 



no output 
requested 



CALL FL0TST3E 

no output 
requested 



computation 

proceeds 



V 



R 


CALL 


PRELIM 


E 

P 


CALL 


SHEETF 


0 

R 


CALL 


SHCONN 


T 


CALL 


ENDNE U 



output requested — CALL NEUDISP 

| 

scans through 
the grid- 
changes coor- 
dinates-find 
node with max 
deformation- 
display first 
card- 



•>— output requested 

i 

CALL NEUSTRE 



I sea ns through 
| the elements- 
| change coord 
I ina tes- 
| com pute Von- 
| Hises stress 
I -display 
ele ment 
results 



I 



7 ' 



RETURN 



scans through 
the grid- 
display nodal 
results 



f 



RETU RN 
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IV. DESCRIPTION OF THE DNDER BATER EXP LOSION 
A. PRESENTATION 

An explosion is a chemical or nuclear reaction in a 
substance (water) that converts an original material into 
other products plus a significant amount of energy. The 
process of tha explosion produces vary high temperatures and 
pressures and occurs with extreme rapidity. As rhe resul- of 
the explosion, the initial mass of explosive becomes a very 
hot mass of gas at tremendous pressures; it is then obvious 
that these conditions will affect the surrounding medium. 

The fact that the water is compressible leads to the 
conclusion that the pressure applied at some location in the 
liquid will propagate through it as a wave disturbance 
[Ref. 9]. It must be pointed out that the pressures 
involved in an underwater explosion are usually so large 
that the wave velocity cannot be assumed independent of 
pressure. Thus, the small amplitude wave theory detailed in 
[Ref. 9] does not apply and this will be the source of many 
complications in describing the behavior of the shock wave. 

The first cause of disturbance to rhe water in an under- 
water explosion is the occurrence of the pressure step at 
the water boundary. Immediately upon its arrival, the pres- 
sure begins to be relieved by an intense pressure wave and 
outward motion of the water. For explosives like TNT or for 
nuclear explosions, the pressure rise can be considered as a 
discontinuous step, and is then followed by a roughly expo- 
nential decay. The duration of the phenomenon is of the 
order of a few milliseconds (more for nuclear explosions) . 
Once initiated, the pressure disturbance is propagated radi- 
ally outward as a compression wave, also called a shock wave 
because of the steep pressure step at its front. 
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Close to the explosion, the velocity of the wave is several 
times the "acoustic" speed of the small amplitude theory 
(c=5000 ft/sec) ; but approaches this limit very rapidly as 
the wave advances outward. 

The pressure level in the wave falls more rapidly with the 
radial distance than what is predicted with the small ampli- 
tude theory, but approaches this behavior at large 
distances. 

B. BOBBLE EFFECT 

As a result of the explosion, the initial mass of explo- 
sives becomes a hot mass of gas at tremendous pressures. 
After the principal part of the shock wave has been emitted, 
the gas pressure is considerably decreased but is still much 
higher than the equilibrium pressure. The water in the imme- 
diate region of the sphere or "bubble" of gas has a large 
outward velocity and the diameter of the bubble increases 
rapidly. The expansion continues and the internal gas pres- 
sure decreases gradually, but the motion persists because of 
the inertia of the outward flowing water. When the gas 
pressure falls below the equilibrium value, the pressure 
defect brings the outward flow to a stop and the boundary of 
the bubble begins to contract at an increasing rate. The 
inward motion continues until the compressibility of the gas 
reverses the motion. Thus, the inertia of the water 
together with the elastic properties of the gas provide the 
necessary conditions for an oscillatory system. The oscilla- 
tions of the gas sphere may persist a number of cycles, ten 
or more oscillations having been detected in favorable 



cases. 



C. SURFACE EFFECTS 

In the case of explosions occurring at shallow depths, 
surface effects will complicate the aforementioned sequence 
of events. When the shock wave hits the surface, the atmos- 
phere cannot supply appreciable resistance by compression. 
As a result, a reflected wave with a negative pressure 
satisfying the zero-pressure condition at the surface is 
formed (figure 4.1). Thus, the resultant pressure observed 
is the sum of the direct and reflected pressures. Therefore, 
the reflected wave arriving at point A will create a sudder, 
drop of the pressure to a smaller value. This is known as 
the M cut-off" phenomenon, typical of free surface effects 
(figure 4.2). 



i 1 




Figure 4.1 Surface Effect on a Shock Have. 



D. PRESSURE DETERMINATION 

As detailed in a previous section, the pressure decay at 
any point is roughly exponential so that it can be written: 

P t (A) = P 0 (A) e^ (4.1) 
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It has been found that the fundamental descriptors cf an 
underwater explosion attack are the charge size (equivalent 
weight of TNT) and the charge standoff (shortest distance 
between charge and target). Theoretical developments about 
the spherical wave detailed in [Ref. 10] have shown that the 
peak pressure P m at any point can be reasonably approximated 
by: 







(4.2) 



where W is the charge size in pouads of TNT and R is the 
standoff distance in feet. 

It has been shown as well: 



V3 V3 a 2 
K 2 H <W /R) 2 



(4.3) 



K l rK 2 ,A l ,A 2 are empirically determined factors that depend 
on the type of explosives used. Iheir values for several 
types of explosives are given on table V. 
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Figure 4.2 Cut-Off Phenomenon. 



TABLE V 

Explosion Parameters 





HBX-1 


TNT 


PENT 


NUKE 


Pnrax 


22347.6 


22505 


24589 


4.32X1C 0 


*i 


1.144 


1.18 


1.194 


1.13 


Decay Constant 


.056 


.053 


.052 


2.274 


A 2 


.247 


.185 


.257 


.22 



E. THE EXP LOS I OH IN EPS! 

EPSA features two different way of describing underwater 
explosions: 

-A discrete form in which the user inputs the pressure 
history at a finite number of times. 

-A functional form that uses equations 4.2 and 4.3 .The 
program requests then the various coefficients and parame- 
ters describing the explosion. 




TIME (MSEC) 



Figure 4.3 Incident Pressure Decay. 
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The presence of free -surface effects can also be 
accounted for with the input of a out-off time after which 
the incident pressure is set to zero (figure 4.3) . 
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V. ANALY SIS AND RESULTS 



A. HODELS STUDIED 

In order to compare the results of the numerical anal- 
ysis with the experimental data, all attention was focused 
on the Explosive Power Meter (EPM> model for which field 
tests had been conducted. The SPM model is a stiffened 
cylinder with end caps whose dimensions are given in figure 
5.1 . 



S'-IO" 




EPS A MODEL (CROSS SECTION) 

END RING 

if DISCRETE 
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if 
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Figure 5.1 Explosive Power Meter. 
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Rather than modelling the end caps with additional 
sheets , it was found to be more efficient to take into 
account the behavior of the end caps using two rigid end 
blocks (figure 5.1). The effects of the explosion on the 
deformation of the end blocks is negligible compared to the 
deformation of points located outside of the end blocks. 
Therefore, it was assumed that the morions of the end blocks 
were pure rigid body displacements. 

In order to gain some insight into rhe influence exerted 
by the stiffeners and the end blocks, a preliminary analysis 
was conducted on a ring stiffened cylinder without end 
blocks as well as on an unstiffened cylinder . 

In addition to the study of the EPM model, the influence 
that the location of the stiffeners had on the deformations 
throughout the shell was evaluated. By performing a compara- 
tive analysis of the deformations, it was intended to mini- 
mize and control the damage caused to the structure. 

The cylinders tested were subjected to an explosion 
occurring at the distance R= 200” from rhe cylinder. The 
location of the explosion was symmetrical with respect to 
the longitudinal and transverse axis of the cylinder (figure 
5.2). A spherical type TNT explosion of intensity W=50 lb 
was selected for the study. It was therefore determined by 
the following parameters (chapter IV) : 
k 1 = 1.18 A 2 = - . 185 

K x = 22505 K 2 = . 058 

The symmetry with respect to the xy and yz plane has 
been taken advantage of by modelling one quarter of the 
model. After a certain amount of sensitivity analysis was 
performed on simple grids, a finite element grid consisting 
of 15 17 nodes and 1440 elements was selected (figure 5.3) . 

For each of the cylinders studied, the time step chosen was 
equal to At = 3. 1 0* 6 s . The explosion process was studied 
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Figure 5.2 Explosion Location. 




Figure 5.3 FEH Discretization 1517 nodes, 1440 element s. 
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over a time interval of 800 time steps, that allowed for a 
shock and after-shock analysis. 

The displacemen ts computed by EPSA consist of a combina- 
tion of rigid body displacements and pure deformations. The 
deformation modes are of significant interest since they may 
induce buckling and even lead to the destruction of the 
structure. As mentioned earlier, the very stiff end blocks 
have pure rigid body displacements. The displacement of 
each node of the end block was subtracted from the displace- 
ment of each node of the corresponding row, giving the 
component of the displacement due to pure deformation. 

For each of the aforementioned cylinders, the deformed 
geometry and the color-filled contour plots of von Mises 
stresses as well as normal displacement were displayed. 
Using identical color assignments, the progressive gross 
evolution of the previous quantities were evaluated so as to 
allow for a comparative analysis of the evolution of phys- 
ical parameters throughout the shell. 

Color-filled contour plots allow for a global representation 
of a quantity distribution and have been found extremely 
valuable in the interpretation of the results. 

For printing and processing reasons, it was not possible 
to include color pictures in this document. Instead, the 
contour plots of the physical quantities under study have 
been included. 

B. ANALYSIS OF RING STIFFENED CYLINDERS WITH END BLOCKS 
1 . 2PM Model 

The contour plots of von Mises stresses at time 
steps 20, 60, 100, 140, 180, 200 are provided on figure A.1 
to A. 6 . As the shock wave hits the structure, it appears 
as if the stresses propagate through the shell and reach 
their maximum value fairly quickly in 50 time steps. After 
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100 time steps, when the structure is fully enveloped by the 
shock wave, the region close to the end block becomes 
heavily stressed (figure A. 4). in addition, soma concentra- 
tion of stresses at the locations of the stiffeners can be 
seen (figure A. 5). At later time staps, the pressure becomes 
decreased and there is a relaxation of stresses. However, 
the region close to the end block as well as some spots 
located around the stiffeners remain heavily stressed, which 
may indicate local buckling (figure A. 6) . 

2 • Contro ll ed D eform at ion s 

In order to obtain a more uniform distribution of 
displacements, the stiffeners have been shifted towards the 
end blocks. The time history evolution of the displacements 
was studied over an interval of 803 time steps for the EPM 
model as well as for the model with shifted stiffeners 
called EPM2 . The contour plots of normal displacements at 
time steps 200, 400, 600, 800 for both models are provided 

on figure A. 7 to A. 14 

The EPM model shows a significant concentration of 
deformations occuring, even at late time steps (figure 
A. 10) , indicating a possibility of buckling. Although unex- 
pected, the fact that the region close to the end blocks 
undergoes the most severe deformations has been confirmed by 
experimental data. A possible explanation to this phenomenon 
is that the inertia of the cross-section of the cylinder is 
relatively uniform along the cylinder, except at the end 
blocks where it jumps to a much higher value. This disconti- 
nuity in the inertia results in concentrations of stresses 
that cause the aforementioned phenomenon. 

The cross-section inertia of the EPM2 model 
increases more smoothly because of the distribution of stif- 
feners along its axis. Thus, the concentration of stresses 
has a lower magnitude and the region close to the end block 
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suffers less damage than in the EPM case. It can also be 
seen that the deformations are more uniformly distributed 
along the axis of the cylinder. Above all, the deformations 
at the late time steps are not as large, indicating that the 
chances of buckling are significantly lower for the EPM2 
model (figure A . 14) . 

Therefore, by performing an optimization of the 
location of the stiffeners, the designer can counteract the 
concentrations of stresses and the buckling phenomena that 
occur in the region close to the end blocks. It is believed 
that controlled deformations can have a very significant 
influence on the survivability of a structure submitted to a 
shock wave. 

C. ANALYSIS DF ONSTIFFENED AND RING-STIFFENED CYLINDER 

It was decided to study the progressive gross responses 
of an unstiffened cylinder as well as that of a ring- 
stiffened cylinder without end-blocks. Both cylinders have 
the same external dimensions as the EPM model. The ring- 
stiffened cylinder is similar to the EPM model except for 
the fact that the end-block has bean replaced by a standard 
stiffener. The evolution of von Mises stresses at time 
steps 40, 80, 100 is provided in figures A. 15 through A. 17 

for the unstiffened cylinder. The evolution of von Mises 
stresses at time steps 40, 80, 100, 150 is provided in 

figures A. 18 through A. 21 for the ring-stiffened cylinder 
without end-blocks. 

For the unstiffened cylinder it is observed that the 
stresses propagate quickly throughout the shell and that 
within a hundred time steps an instability phenomenon occurs 
showing the existence of local buckling (figure A. 17). At 
later time steps the buckling spreads over the entire stiff- 
ened shell. 
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The evolution of von Mises stresses in figures A. 18 
through A. 21 shows that the ring-stiffened cylinder can 
withstand much higher stress levels than the unstiffened 
shall. This result was expected since the stiffeners provide 
the stiffness required to withstand higher loads. At time 
step 100 the unstiffened shell is already subjected to local 
instability characterized by a "hard spot" in the middle of 
the model (figure A. 17). On the other hand, the stresses in 
the stiffened cylinder are much more evenly distributed 
throughout the model, with high amounts of stresses concen- 
trated around the locations of the stiffeners (figure A. 21). 

It can also be observed that a significant concentration 
of stresses occurs at the extremities of the stiffened shell 
(figure A. 21). Recalling that the end-block has been 
replaced by a standard stiffener, the cross-section of the 
shell has a greater inertia at the extremities due to the 
facr that the two stiffeners located at the extremities are 
close to each other. Therefore this phenomenon is similar to 
the one encountered when studying the EPM model. However 
the concentration of stresses for the ring-stiffened cylin- 
drical shell is less significant than for the EPM model, due 
to a smaller discontinuity in cross- section inertia. 
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VI. CONCLUSION 



A FORTRAN module that merges the finite element coda 
EPSA with the color graphics system PATRAN-G has been 
designed and succasfully completed. The non-linear elasto- 
plastic responses of various types of submerged cylindrical 
shells have been evaluated using the EPSA/PATRAN-G system. 

The analysis of the progressive gross responses of a 
ring-stiffened cylindrical shell with end-blocks is believed 
to provide useful information regarding the behavior of a 
submerged structure subjected to an underwater explosion. 
The influence of the location of the stiffeners on the 
deformations has been studied and is also believed to be of 
significant help in the determination of an optimal design 
that will minimize the damage due to an underwater 
explosion. 

The utilization of the color graphics system in the 
interpretation of the results of analysis has been found to 
be an extremely valuable tool. It is the author's belief 
that the use of color graphics systems will become increas- 
ingly important in the analysis of complex phenomena such as 
underwater explosions on submerged structures. 
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AP PENDIX A 
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Figure A. 1 



Von Bises stresses. 
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Figure A. 2 Von Mises stresses, time step 60 
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Figure A.3 Ton Mises stresses, time step 100, EPH model. 
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Figure A. 4 Ton Mises stresses, time step 140, EPfl model. 
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Figure a. 5 Von Bises stresses, time step 180, EPM model. 
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Figure &. 6 Von His es stresses, time step 200, EPH model. 
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Pigure A. 7 Horaal Displacements, step 200, EPH model. 
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Figure A. 8 Noraal Displacements, step 400, EPH model. 
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Figure i.9 Hormal Displacements, step 600, EPS model. 
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Figure A. 10 Normal Displacements, step 800, EPH model. 
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.11 Hormal Displacements, step 200, EP M2 model. 
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Figure A. 12 Normal Displacements, step 400, EPH2 model. 
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Figure A. 13 Horaal Displacements, step 600, EP M2 model. 
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Figure A. 14 Hormal Displacements, step 800, EPH2 model. 
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Von aises stresses, step 40 , unstif. shell. 
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Figure 4.15 



UNSTIFFENED CYLINDRICAL SHELL 
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Figure A. 17 Von flises stresses, step 100, unstif. shell. 
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Figure A. 18 Voa Hises stresses, step 40 , stiff, shell. 
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Figure A. 19 You Hises stresses, step 80 , stiff, shell. 
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Figure A. 20 7on Mises stresses, step 100, stiff, shell. 
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Pigure A. 21 Von Hises stresses, step 150, stiff, shell. 
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APPENDIX B 

REVIEW OF NONLINEAR FINITE ELEHENTS 

A. IHTRGDOCTION 

This appendix is intended to give the reader seme 
insights into nonlinear finite elements. The reader is 
assumed to have some previous knowledge of finite element 
theory. The basic principles of the theory will be quickly 
reviewed, but the study will focus on the problems that 
occur when dealing with nonlinear theory. Most of the 
information has been taken from [Ref. 6 ] as well as from 
the course the author had at M.I.T. with K.J. 3athe in 1982. 

B. THE HEED FOR A NEW THEORY 

Considering a coordinate frame defined by (i,j,k), a 
body of volume V in which point A(xi,x 2 ,x 3 ) is subjected to 
the displacements (U]_,u 2 ,u 3 ), corresponding to a strain 
vector {e} (figure B.1). 

In the following sections, the superscripts 0 and t will 
refer to the body at time 0 and t respectively, the 
subscripts 0 and t will refer to the configuration at time 0 
and t respectively. This chapter, for the purpose of 
simplicity, will first deal with the static nonlinear anal- 
ysis of the material. 

In the linear theory of finite elements, one uses the 

well known Cauchy stress tensor associated with the 

engineering strain tensor e^ =_ff3Ui + 3u^l . 

2|_3 x j 5ciJ 

Then, the principle of virtual work is written: 

where represents the virtual work of externally applied 

forces and Se is a virtual strain. 
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Figure B.1 Geometric Conventions. 



i tSe mn = (B.1) 

mn 11111 

Equation (B.1) is then discretized over the body and becomes 
a set of integrations over each of the finite elements. In 
the case of a large displacement, the volume of the body 
over which the integration is performed might have signifi- 
cantly changed. Also notice that equation (B.1) is written 
in the original coordinate frame (defined at t = 3) and that 
and refer to the current configuration of the body. 

The Cauchy stresses at time t+ a t cannot be obtained by 
adding an increment due to the straining of the material to 
the stresses at time t . The rigid body rotation of the 
material has to be taken into amount since the Cauchy 
stresses vary under rigid body motions. Therefore, we must 
perform the integration of equation (B. 1) over the unknown 
current volume with respect to the original geometry that 
could be significantly different from the current one. 

The above discussion emphasizes the need for a new set 
of stress and strain tensors that would alleviate the afore- 
mentioned problems and enable the integration of the prin- 
ciple of virtual work to be performed. 
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C. DEFINING NEW STRESS AND STRAIN TENSORS 



1 . Green- L a gr ange Str a in T enso r 

The structure introduced earlier is considered and 
two points A and B at t=0, of coordinates^ (°x) , °B(°x) 

are defined . At time t r the body has been deformed and A 
and B have moved to t A( t x) and t B( t x) (figure B.2). 

C { 




Figure B.2 Displacements Conventions. 

A Taylor expansion is used to express the coordinates of B 
as a function of the coordinates of A. 



\ \ ( 0 Xj- ^ 

3 



(B.2) 



or with d t x i = t x i - t x i and d °x L = °x i - \ , gives : 



(d V) = ( 3_^ i ) (3 V) 



3 X n 



(B .3 ) 



{ d fc x} = [ ] {d°x} 



(3.4) 



where [ &C ] = ( ) , (<14} = (d fc xi) , {d °x} = (d°Xi) 

In other words, aquation (B.4) expresses how a small fiber 
defined by the vector {d °x} at r = 0 has rotated and 
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extended between time 0 and time t when it becomes {d fc x} . 
The matrix [gX] is called the "deformation gradient" . 

The new length dS of the fiber will be : 

( t dS) 2 = (d t x} T {d fc x} ( T refers to the transposed 
matrix) 

and therefore 

<^S)2 = (C^X] [d 0 X}) T ([gl] {d°x} ) 



( fc dS)2 = {d °x } Co c 3 {5°x} 



(B.5) 



t t T t 

[ 0 C ] = C o x 3 [ o x ] is called the "Cauchy-Green deformation 
tensor" . 

Notice that if the Cauchy-Green deformation tensor is iden- 
tity, equation (B.5) indicates that the length of any fiber 
will not vary. In other words, whenever rigid body motions 
are considered, the Cauchy Green deformation tensor is iden- 
tity since the fibers do not stretch. 

The principles of the finite element method will now 
be quickly recalled. Assume that the displacement ) of 
any point of the body can be written as a combination of the 
displacements of N selected points called "nodes" : 



<v * £ w' 

k-l 



(B.6) 



Where the h are interpolation functions that depend only on 
the geometry of the body. In addition, the nodes are chosen 
so as to get a division into quadrilateral elements and it 
is assumed that the displacement of any point is only a 
function of the displacements of the corner nodes of the 
element it beLongs to. Then equation (B.6) becomes 

4 



( 



Ui) = £ K {U i 



) J 



(B .7) 



k=l 



Recalling that (u^ is simply ( fc xi) - (°xi) gives: 
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(B.8) 



4 



( c xl = ( X) 




k=l 



The deformation matrix [ X] 
using the previous equation : 



matrix X] can be expressed very simply 




(B.9) 



Define now the "Green Lagrange" strain tensor by : 



3 * V2 «$C] - [I]) 



(B. 10) 



Where [I] is the identity matrix. 

From the previous derivations, it can be observed that the 
Green Lagrange strain tensor is 0 for rigid body morions. 
The Green-Lagrange strain tensor refers to the body at time 
t with respect to the initial configuration. This is why it 
will be so useful in dealing with large deformations. 



ciple of virtual work to the structure under study. In 
particular, having defined a new strain tensor, the relation 
giving the virtual Green-Lagrange strain tensor corre- 
sponding to a virtual displacement (6 u^ must be known. In 
the case of a linear problem, the virtual engineering strain 
tensor would be: 

fit 2 ;- = 1/2 ( 3ui + 3Uj ) 

3 3 x j 3 jg. 

In the case of large displacements, the Green-Lagrange 
strain tensor should be used. It is shown in [ Bef . 6] that 
the virtual Green Lagrange strain tensor corresponding to 
virtual engineering strains is: 



Recall that the ultimate goal is to apply the prin- 





(B. 11) 



or: 
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(B. 12) 



$ t®i j = ^2Li ' J x j ^0 £ i j 
3x m 3x n 

Having defined a strain tensor which is invariant under 
rigid body motions, a stress tensor corresponding to the 
Green-Lagrange strain tensor needs to be defined. 

2 . Str ess Mea su res 

Starting with the Cauchy stress tensor , the 

Piola-Kirchoff stress tensor is defined : 

O^ij t^jjn tT mnt x jn (B.13) 

Where °p , are the densities of the material at time 0 and 
t respectively and ( t x i^ = [ tX 3 1 is the inverse of the 
deformation tensor defined previously. 

Equation B. 13 can be easily rewritten in equivalent form: 

"'mn = -^ O x in n (B.14) 

It can also be shown by using the principle of mass conser- 
vation that : 

Op = t p det[5x] (B. 15) 



D. PRINCIPLE OF VIRTUAL WORK 



The principle of virtual work of equation (B.1) is 
rewritten using the strain and stress tensors defined previ- 
ously in equations (B.12) and (B.14) : 




0 ^Lk O^k § x im ? x j§ 0 £ ij ^ 



(B. 16) 
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(B. 17) 



and 




using equation (B.15) : 






oSij6oei jdV 



(B. 18) 



Thus, the principle of virtual 
expressed in terms of a new sat 
tensors, integrated over the origina 



work has 
of stress 
1 volume °7 



been 

and 



simply 

strain 



E. THE IHCRESEHTAL CONTINUUM MECHANICS EQUATIONS 



In this section, the principle of virtual work will be 
applied to tha structure and tha incremental formulation 
using the Piola- Kirchcf f and Green- Lagrange tensors will be 
developed . Non linear terms will arise from the rather 
complicated definition of the strain tensor, but it must be 
pointed out that the new formulation provide the means for 
the modelling of large deformations. 

Assume that the configuration of the body at time t is 
known, the configuration at time t+At must be determined. 
Writing the principle of virtual work at time t+At gives : 




t+ ^R : external work at time t+At 
6 o e ij virtual increment in G.L. strain 
o^ij : stress state at time t+At 



Separating between the known tarms that refer to the config- 
uration at time t and the unknown terms that are the incre- 
ments of stress and strain between time t and time t+At 
gives : 
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t+ At t t 

0 s = o s ij + oSij and e = o £ ij + 0 £ ij e (3.20) 

oSij ando £ ij are simply the increments in stress and strain 
respectively. 

It is derived in [Ref. 6] that the increment in 
Green-Lagrange strain o £ ij is made of a linear part oeijand a 
nonlinear oneji^ . The term linear refers to the increment 
in displacement u^ . 



o e ij = c%j + o n ij cr 6 = 6 ( oeg + onij ) (B.21) 



Using equation (B.21) in equation (B. 19) gives : 
[ ($5 y + oSij ( 5 (oeij + 0 0ij )) iV = ^ 

•V 



(B. 22) 



Again, separating between the known and unknown terms gives 



1° 



^ - o n ij 



t+< 






jt 'U 5Lr 



oeij 



(3. 23) 



For small increments in displacements, equation (B.23) 
written at time t indicates that 5 0 eij = 6o£ij • This signi- 
fies that in equation (B.21) the non linear term is negli- 
gible 

Then the constitutive law of the material detailed in 
chapter II allows to relate stresses and strains: 



O^ig “ O^ijrsC^ij ^ 0^-ijrs 0 e i 



(B. 24) 



and B.23 becomes: 



[& o«Lj cfcjsdfad* * [ SSg6 0 r^dV= - ( <Aj Cfij* v 

J o v *°V -°v 



(B. 25) 
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The right hand side terms of the previous equation are 
known : 

is the linear increment in virtual strain involving only 
known terms 

Bsij is known from the previous time step 

t>At . 

R is the work of external prescribed virtual forces. 

The left hand side of B. 25 is unknown since 0 e^ja nd 
involve the displacement from time t to time t+At . 

F. FINITE ELEMENT DISCRETIZATION 

Equation (B. 25) will be discretized over the structure, 
using the finite element approximation defined in section 
B.9 . Let N be the number of nodes, the principle of 
virtual work will be invoked N times, setting a unit 
displacement at each node in turn: 

o u k = 1 , <$Uj=0 , k*j 

A system of equations whose unknowns are the nodal displace- 
ments ' is obtained Let { A u} be the vector of unknown 
displacements, {F} be the vector of nodal point forces 
equivalent to the internal stresses. 

Then equation (B.25) can be rewritten in matrix form : 

[Jk d ]{Au} + [% L ]{Au} = (B. 26 ) 



CqK l ] and t qK^] are known from the material char acterist ics 
at time t and correspond respectively to a linear and nonli- 
near contribution. It is therefore possible to solve equa- 
tion (B.25) for {Au} . However, because of the assumption in 
equation (B.24), the exact solution might not be reached 
immediately. Furthermore, depending on the time step size, 
the solution process might even be unstable! In any case. 
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an iteration for solving aquation (B.26) musz be performed 
until the exact solution is reached . 

A widely used scheme is the "modified Newton iteration" 
defined by the following equation and boundary conditions : 



( [Sk l ] + [JknlJ ) (Au}i 



= f 1 



(B. 27) 



and 1 f^ t u} i-; ^ {Au]^ • with the initial conditions : 

n 4- teA: o t 

l u}° = (M and { F}° = C?} 

£AU} 1 is the vector of incremental nodal point displacements 

at iteration i. 
ta-At 

{ R} is the vector of applied loads (constant in the itera- 
tion) 

*t£ A 

{ F) is the vector of nodal point forces equivalent to the 
stresses at time t+A t, iteration i-1 . 

At the first iteration, equation (3.27) reduces to equa- 
tion (B.26) giving an increment of displacement {Au} . Then 
a better approximation of 0 e rs • Scfii i s obtained . gSri i s 
updated to the new state of stresses and becomes o^ij 
Equation (B.27) is then used to determine the new increment 
in displacement (Au) 
displacement is 
equation (B.27) 



, and so on until the increment in 

t3- fc. tHAt i 

displacement is small enough, so that { 3} = { F } in 



G. INCLUSION OF DYNAMIC PORCES 

If the loads are applied rapidly, inertia forces need to 
be considered and a truly dynamic problem has to be solved. 
Using d’Alembert's principle, the eLement inertia forces are 
simply included as part of the body forces. Let {u} be the 
vector of nodal accelerations and [ M ] be the mass matrix of 
the system. Then the principle of virtual work is written in 
the following way : 



85 



R -CM} (au) 



ta-Zt 
- F 



(B . 28) 



IcSyuu) * tS^tau ) 



t3-& 



Equation (3.28) represents a system of differential equa- 
tions of second order. If the non linear term (^k KL ] were 
negligible, the solution could be obtained by standard 
procedures for solving differential equations with constant 
coefficients. However, the procedures proposed for the solu- 
tion of general systems of differential equations can become 
very expensive if the order of the matrices is large. 
Therefore, whenever the system is linear or nonlinear, seme 
effective methods for solving the equations governing the 
equilibrium are required. 

1 a Direct Integ ration Metho ds 

The essence of direct integration methods is based 
on two ideas. First, it is aimed to satisfy B. 28 only at 
certain time intervals apart. Second, a variation of accel- 
eration velocities and displacements is assumed within each 
time step. The form of the assumption determines the accu- 
racy, stability and cost of each method. 

In the following, assume that the initial conditions 
(displacements, accelerations, velocities) at time 0, 
denoted ( °u ,°u , u ) are known. In the solution, the time 
span under consideration, T , is subdivided into n equal 
time intervals At. Assuming that the solution is known at 
time t , the methods of getting the solution at time t+At 
will be investigated. 

2 . Centr al Diff erence Met hod 

In the Central Difference method, a finite differ- 
ence approximation will give the acceleration at time t : 



t •• i 

U — "r 

At 



t- At 

2 ( u 



-,t tiAt 4 

- 2 u + u ) 



(B. 29) 
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Writing the principle of virtual work at time t and 
substituting into equation (B.30) gives then: 



m^u} + CoK ]{ fc u} = { fc R} 



(B. 30) 




= { fc R} 



t ? t , tfit 

(Co K] - 2[»]/At 2 ) { u} ^ra]{ u} (B . 31 ) 



The previous equation gives the deformation at time t+At 
from the characteristics of the system at time t . 

When [ M ] is diagonal, which is frequently the case 
for mass matrices, the solution at time t+ At does not 
involve any triangular factorization of the matrix [ M ] , 

thus leading to more efficient computations. 

The shortcoming in the use of the central difference 
method lies in the time step restriction: for stability, the 
time step size t must be smaller than a critical time step 
/$ cr which is equal to T , where T is the smallest period 
cf the finite element system. 

The central difference scheme is fairly easy to 
implement for the integration of a system of non linear 
differential equations. However, because of the limitations 
of the time step , it might not be suitable for cases when 
loads are varying at a slow pace. 

3 • Impli cit Int eoratio n Sche mes 

Since the interest of this study lies in central 
difference schemes, this section will be limited to a short 
description of the fundamentals of implicit time integration 
schemes. 
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Implicit time integration schemes use the principle 
of virtual work written at time t+ At and not at time z as 
for the central difference method : 

tJ-At t-tAt tnAt tg-At 

[ M] {u} ♦ CqK] {u} = { R} (3.32) 

taAt 

Again, using a finite difference approximation of {u} and 

. . . tgAt 

replacing into equation (3.32) enables to solve for { u} . 

fetAt 

Since the formulation involves the rigidity matrix fo K ] and 

■tf-At 

the external work { R} which are both unknown, the system 
has to be solved in a similar way to the Modified Newton 
iteration that was detailed previously. 

Implicit time integration schemes are stable regard- 
less of the size of the time step used. However, if the time 
step size is too large, significant errors can be accumu- 
lated at each time step, leading to unrealistic results. 

The reader will find more details on the various 
implicit method in [Ref. 6], let, it can pointed out that 

implicit methods are more tedious to implement. On the other 
hand, a larger time step can be used in the solution proce- 
dure, which can be of extreme importance when studying 
phenomena over a significant period of time. 
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APPENDIX C 

HOW TO USE THE THANSLATOR EPS A-PATR AN-G 

This appendix is intended to explain in derail the use 
of EPSA/PATRAN-G post-processing facilities. It is divided 
in three sections that will deal successively with 1) the 
displaying of the original model ; 2) the deformed geometry; 

3) the contour plots of element and nodal points quantities. 

A. DISPLAYING THE ORIGINAL MODEL 

When making an initial EPSA analysis on a particular 
structure, the geometry of the model has to ba input into 
PATRAN-G. As explained in chapter III, all the geometrical 
information is contained in a file FOR019.DAT that is 
created each time an EPSA run is made. The input of the 
original geometry must be made via the neutral input mode of 
EPSA. The procedure, starting from the "logon" to PATRAN-G 
is the following ; 

- Select the 30 option 

- Select the new data file option (option 1) 

- Select the neutral system (option 4) 

- Select the input mode (option 4) 

- Input the neutal file name : FOR019.DAT 

Tha original geometry will then be displayed on the screen 
It is often found convenient to have a perpective view 
of the model under study. In that case, the user should : 

- Issue the VIEW command 

- Select the rotation about the absolute axes (option 1) 

- Input an angle of rotation 

(23,-34,0 will give a very nice view but any angle can be 
input) 
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- Issue a PLOT command to have the model displayed in the 
new a xe s . 

An example of the procedure is provided on table VI. 



: 



TABLE VI 

Finite Element Model Input Procedure 



MODE? 1. GEOMETRY MODEL 2. ANALYSIS MODEL 3. DISPLAY 4 . NEUTRAL SYS. 5. END 
>4 

NEUTRAL FILE? 1. CREATE OUTPUT 2. INPUT MODEL 3. POST-PROCESSING 4. END 
>2 

INPUT NEUTRAL FILE NAME 
✓FOR019.DAT 

DO YOU UISH TO OFFSET ANY NEUTRAL INPUT IDS? (Y/N) 

>N 

EPM 200 STEPS , NO STIFFENERS U*30. 

SHALL UE PROCEED UITH THE READING OF THIS FILE 7 (Y/N) 

> 






The PLOT command can be issued anytime to display the orig- 
inal geometry on the screen. 

When studying complex models, one does not want the 
element and node numbers to be printed along with the geom- 
etry of the structure. The command SET, LABS 3, OFF followed 
by a PLOT will display the original geometry without any 
labels printed. 

When a model has been input into PATRAN-G, a data file 
PATRAN.DAT is created on the user's directory. When 
connecting with PATRAN-G at a later time, the user can 



90 



select the option "last data file " (option3) to have the 
original geometry displayed on the screen without having to 
input the model again. 

B. DEFORMED GEOMETRY 

On the second card of the PATRAN-G input deck, the user 
specifies the number of PATRAN-G displacement data file and 
the number of element results data file. Two non-zero inte- 
gers in free format must be placed at the end of the second 
card (section II in the user's manual) requesting the number 
of displacements and of elements files respectively. 

Assuming that the user has made a 200 steps run with 10 
output requests for PATRAN-G displacement files, ten (10) 
files FOR018.DAT will then be created at equal time inter- 
vals. The deformed geometry corresponding to time step 100 
will therefore be contained in FOR01 8 . DAT; 5 . 

To display the deformed geometry corresponding to time 
step 100, the user should issue the following commands: 

- RUN,DEF : requests deformed geometry option 

- Input the name of the displacements file : FORD 1 8 . DAT 5 

- Select the PLOT option (option 3) 

Select the undefcrmed geometry (2) followed by the 

deformed geometry (3). An example of the procedure is 

provided on table VII. The undeformed geometry superposed 
with the deformed geometry will taen be displayed on the 
screen. 

C. POST-PROCESS IHG OF ANALYIS RESULTS. 

Element-related quantities like von Mises stresses are 
contained in POR016.DAT files, nodal point quantities are 
stored in FOR018.DAT files. As described in chapter III, 
each column of those files contains a specific quantity 
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TABLE VII 

Deformed Geometry Procedure 

> 

MODE? 1. GEOMETRY MODEL 2. ANALYSIS MODEL 3. DISPLAY 4. NEUTRAL SYS. 5. END 
>RUN, CON, COL, 4 

CURRENT FILE FOR NODAL RESULTS IS FOR018.DAT 
NEUTRAL RESULTS FILE? l.NEU FILE 2. CURRENT FILE 
>1 

INPUT THE RESULTS FILE NAME: 

>FOR018.DAT;5 
DATA UIDTH * 5 

FILE TITLE »EPM 800 STEPS 

7.5000129E-04 



DATA UALUES RANGE FROM -0.525E+00 TO 0.134E+00 
ASSIGNMENT? l.AUTO 2. MANUAL 3. SEMI-AUTO 4. USE CURRENT LEUELS 5. END 
>4 

PREUIOUS CONTOUR LEUELS USED. 

MODE? 1. GEOMETRY MODEL 2. ANALYSIS MODEL 3. DISPLAY 4. NEUTRAL SYS. 5. END 
>RUN,HI,C 



(i.e. column 31 of FOR016.DAT contains the von Mises 
stre-sses). The reader will refer to chapter III for the 
detailed organization of those files. 

Again, assume a 200 time steps analysis, with 10 output 
requests for element results files. The user might want to 
display the contour plots of von Mises stresses at time step 
100 . In this case the following commands should be input : 

- RON, CON, COL, 31: tells PATRAN-G to look at column 31 

that contains the von Mises stresses. 

- Input the file name FOR016.DAT 5 

- PATRAN-G will then ask for a color assignment (automatic, 
manual, semi-automatic, current levels used) that the user 
will select according to his needs. 

The contour plots are then ready to be displayed : 

- RUN, HI, FR would display the color-filled contour plots 

- RUN, HI, CON would display the contour plots (color lines) 
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The von Mises stresses are the most useful element quan- 
tities to be displayed, but other element-related quantities 
detailed in chapter III like the x and y stresses in local 
or global coordinates could be displayed as well by looking 
at their corresponding column in the element results data 
file. 



Dealing with nodal point quantities, the displacement 
normal and the velocity normal to the shell are very mean- 
ingful quantities in an analysis. They are respectively 
stored in column 4 and 5 of the displacement results files 
FOR018.DAT . 

A contour plot of the normal displacement at time step 100 
would then be obtained via the following commands: 

- RON, CON, COL, 4 (look at column 4) 

- FOR018.DAT 5 (name of the file) 

- Color assignment chcsen 

- RON, HI, C or RON, HI, FR 

Not ice When the fluid-structure interaction is ON, the 
normal displacement contained in column 4 corresponds to 
pure defor m ati ons, the rigid body contribution having been 
taken out. 

All the element results processing is implemented in the 
routine NEOSTRE, all the nodal points processing is imple- 
mented in routine NEODISP. It should be pointed out that any 
modification to the capabilities of the translator (i.e. 
being able to display other types of quantities) can be 
mads by modifying these routines only. 
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APPENDIX d 
LISTINGS 



This appendix contains the various files than constitute 
the translator module. The submodule that displays the orig- 
inal geometry is listed on the first four pages. It is 
imbedded in file FRANK. FOR. The submodule that takes care of 
the post-processing facilities is imbedded in file DISP. FOR 
and is listed in the remaining pages. It has been mentioned 
previously that EPSA had been slightly modified to accomo- 
date color graphics capabilities. The only interaction of 
EPSA with the translator occurs via subroutine calls. All 
the "calls" occur in COMPOTE (for post-processing) and 
REPORT (for the original geometry) . A labelled COMMON 
called FRANK has been created and is defined in the routine 
AAA as reguired by EPSA. The requests for PATRAN-G outputs 
are echoed in the EPSA output file, all the modifications 
for that purpose having been made in routine READIN. The 
user must be aware that the size of blank COMMON array A has 
been increased to store the deflections in x, y, z direc- 
tions instead of only the z direction previously. 
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Subroutine dppI im 
dimension A(l) 

Common i a ( 1 ) 

eauivalence (ia(l)#a(l)) 

conmon/ssi ze/ i bq ( 1 ) , nq j , ne 1 tot»nlbd#nload»nbpect/ 

1 nbouad, isheet, norors^ nsots, nstrotS> nvots, nhpts, 

2 nsstyo, nn j , nntot # lqdso, lioud, Ibcalc 

C common /cDara/ nsteos, nbeg, nend» nsheet 

common /stab / ibt(l),jssi2e,jsDar,jvelo,jstre,jxmas,jielm,jbmat, 

1 j 1 1 Dd , j 1odo,joret,j1his, tsfrn, i f o r c , J x l oc , j nq i , j n n i , j nab eg , 

2 j 1 s i de 
c 

common/*rank/nfntot, 1 1 u 

COM*on/titrE/ntitlE(80) 

ch a r a c t er *9 buff 
character's title 

CHARACTER*8 TIM 
CHARACTER *9 VER 
c 

C this routine comoutes new number of nodal points 

c for any given sheet 

C 

KJsJNQBEG-JNNI 

nf n tot =nntot-I A ( JNNI) -I A ( JNQ8EG-1 )-2*(KJ- 2) 

1 1 t ype=25 
1 1 k c = 1 

II i v = 0 

1 I id = 0 
1 1 n 1 =0 
I ln2 = 0 
1 1 n3 = Q 
I 1 n4 = 0 

I I n5 = 0 

w r i t e ( I 1 u , 1 0 ) 1 I t y 0 e , 1 1 i d , 1 1 i v , 1 1 k c , 1 I n 1 , 1 1 n 2 , 1 1 n 3 ^ 1 1 n 4 , 1 ln5 

1 0 format ( i 2/ 9i 8 ) 
write(llu,ll) (NTITLE(I) , 1=1,80) 

11 format(80 Al) 
c 

c talcing care of second oacket 

c 

1 1 t ype = 26 
1 1 kc = 1 
1 lnl=nfntot 
I I n2 = ne 1 t o t 

wri teH 1 u , 1 0 ) I 1 t y o e , 1 1 i d , 1 1 i v , 1 1 k c , 1 Ini, 1 ln2,l 1 n 3 , 1 1 n 4 , 1 ln5 
call da t e ( 8UFF ) 
write(llu,l2) BUFF 

12 format ( A,3X, • 17: 12:09* ,5X, * 1 .4 ' ) 
return 

end 



subroutine shconn 
dimension A(l) 
commo n i a ( 1 ) 

eauivalence Cia(l),a(l)) 

common/ssi ze / ibg(l),naj ,nel tot,nlbd,n1oad,norect, 

1 noauad, isheet, norots, nsots, nstrots, nvpts, nhpts, 

2 nsstvo, nn } , nntot, I gdso , lioud, Ibcalc 
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c common /coara/ nsteo?» nbei/ nendr nsHeet 

common /stab / ibtfl)*jssize*jsoar#jve1otjstretjxmasfjie1mtjbmatr 

1 jllcdrjlodD^joretr i1hisrjstrn,iforcrixloc#)nqi/jnni,jnaDeq/ 

2 j 1 s i de 
c 

ComTton/f rank/nfntoN 1 lu 
c 

c this routine will write the element corner nodes of each 

c element on neutral file 

c 

1 1 t ype = 2 
1 1 kc = 2 

1 1 nod = 4 
1 1 i v=4 
LLN1 =0 
L L N 2 = 0 
LLN3=0 
LLNUrO 
LLN5=0 
LLCONFcO 
LLPI0=0 
LLC£I=0 
THET1=0. 

T HE T2-0 . 

THE T3 = 0 . 

INITIALIZATION DONE 

nor ev = 0 
1 1 e 1 = j nn i - j no i 
c 

do 20 0 k = l, 1 lei 
1 1 row = lA(]nai * k - 1 ) 
c number of elts in each row 
do 10 0 j = 1 , 1 1 row 
1 el - I A ( j nqbeg + K- 1 ) + j - 1 
nlel l=j+norev 
n1el2-i+l+nprev 
nlel 3 = j tltnorev+1 1 rowM 
n1e14=j>norev*1 1 row+1 
c ready to disolav oacket 
c 

1 1 i d-LEL 

writeCl 1 u # 8 0 ) 11 tyoe, 1 1 id, 1 1 i v, 1 1 kc , LLN1 , LLN2 , LLN3 , LLN4 , LLN5 

80 f o r ma t ( i 2 , 8 i 8 ) 

writedlurBl) 11 nod , LLCONF , LL* 1 0 , LLCE I , THE T 1 * THE T 2 , THE T 3 

81 f o rma t ( i 8 , 3 i 8 , 3e 1 6 . 9 ) 

writeCl 1 u , 82 ) nlel I,n1e12,n1e13,n1e14 

82 f o rma t ( 1 0 i 8 ) 

100 continue 

c 

norev = norev+1 1 row+1 
200 continue 
r e t u r n 
end 



sub rout i ne shee t f 
dimension A(l) 
common i a ( 1 ) 

equivalence (ia(l),a(l)) 
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o o r> 



common/ssi z e / i bq ( 1 ) * no j * ne 1 tot»nlb 1 »nloa 3 rnbr*ct» 

1 noauadr isheet/ norotSf nsots# nstPDts» nvots* nnots* 

2 nsstyo* nn j , nntot* I ;dso^ liqud* lbcalc 

c common /co ara/ nsteos* nbeq* nend* nsheet 

common /stab / ibt(l)*jssize*jsoar,jvelo*jstre*jxmas*jielm*jbmat* 

1 j 1 1 od * J lodOr joret f j Ihisr Jst pn f j forcM’ ocr jnoi,jnni,jnabeq* 

2 j 1 s i de 
c 

Comfnon/f pank/nfntot * 1 1 u 
character*! qtyoe 
c 

c renumbers the nodes* chanae coordinates* d i s o 1 a y oacket l 

1 1 t yoe= 1 
1 1 i d=0 

1 1 i v = 0 
1 1 kc = 2 
llnl=0 
1 ln2 = 0 
! 1 n 3 = 0 

i l n a = o 

1 1n5 = 0 
ncoun = 0 
1 1 do f = 6 
1 1icf = l 
gt y oe= ' G ' 

1 1 c on f = 0 
! 1 c i d = 0 
LSPCI=0 
LSPC2=0 
L3PC 3 = 0 
LSPC4=0 
LSPC5=0 
LSPC6=0 

INITIALIZATION DONE 
k j = j nqbeq- j nn i 

k row = k j -2 

do 20 0 j = l * k row 
C looo on new nb* of row 

c 

ncDun=ncoun + IA( jnni + j - 1 ) 

1 1 o t = I A ( jnni ♦ j - l ) - 2 
do 100 1=1 * 1 lot 

11 i dsn id + t 

x l loc = A(ncoun*2 + jxloc+’2*l ) 
yl 1oc = A(ncoun*2tj x 1 oc + 2* 1 + 1 ) 

C skio first node of row 

c 

2 l 1 OC =0 . 

if (A(jvelo-2).ne.O.) then 
rcouro=l ./A(jvelo-2) 
theta=xlloc/rcouro 
*11 oc r rcourb*sin( theta) 

2 lloc=rcourb*cos(theta)-rcourb 
end i f 
c 

if ( a ( j ve 1 o* l ) . ne . 0 . ) then 
rcouro=l . / a ( jvelo-l ) 
thet a=v l loc/rcouro 
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c 



70 



c 

c 

c 



80 

81 



82 

100 

200 



c 

c 



80 



yl 1 oc=rcouro*si n(theta) 
z 1 1oc=rcouro*cos(theta)-rcourb 
end i f 

if((a(Jvelo-l).ne.0.).and.(a(jve1o-2) .ne.O.) ) then 
write(llu#70) 

format (' error, two curvatures are non zero’) 
endi f 

ready to disolay oacket 

wri t e( 1 1 u,80) 11 tyoe, 1 1 i d , 1 1 i v , 1 Ike , 1 1 n 1 , 1 1 n2 , 11 n 3 , 1 1 n« , 1 1 n5 
f o rma t ( i 2 * 8 i 8 ) 

wr i te ( 1 1 u; 8 1 ) 1 1 1 oc t v 1 1 oc , z 1 1 oc 
format ( 3e 1 6 ♦ 0 ) 

wri te(l 1 u r 82 ) 1 1 i c f # qt yoe # 1 1 do f * 1 1 con f # 1 lcidfLSPCl#LSPC2# 
1 LSPC3 LSPCa,LSPC5,LSPC6 

format(Il#alf3i8,2x/6i 1) 
conti nue 
conti nue 
return 
end 



subroutine enaneu 
Common/franK/nfntot/1 lu 
1 1 t y pe=99 
1 1 id=0 
1 1 i v=0 
11 kc = l 

writeCl 1 u # 8 0 ) 1 1 t y o e * 1 1 i d # 1 1 i v # 1 lkc 

format ( i 2 1 3 i 8 ) 
r e t u r n 
end 
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REPORT 
BLANC 2 
BLANC 3 
8LANC 4 

CPARA 2 
CPARA 3 
MSPARA 3 
SSIZE 2 
SSIZE 3 
SSIZE 4 
REPORT 8 
STAB 2 
STAB 3 
STAB 4 
STAB 5 
STAB 6 



INTEGER NSUBT1 (80) ,NSU8T2(80) 

KJ=JNGBEG-JNNI 

NFNT0T = NNT0T-IA (JNNI) -I A ( JNQBEG-1 ) -2MKJ-2) 

LLU = 1 R 
LLN= 1 8 
LLM=17 
DEFMAXsO. 

MAXNO0=NFNT0T 

NWIDTH=5 

NOPAX =0 

LL I 0 = 0 

NCOUN=0 

KROWsKJ-P 

00 10 J= 1 / 80 
NSUBT 1 ( J) =0 
10 NSU8T2(J)=0 

PUT THE TIME STEP ON FIRST SUBTITLE: 

CLOSE(UNIT=LLM) 

OPEN(UNIT=LLM,STATUS='NEW' ) 

HRITECLLM,*) TIME 
CLOSE (UNIT =LLM ) 

RE AO ( LLM , 99 ) (NSUBT 1 ( J ) , J = 1 ,80) 

CLOSE(UNIT=LLM) 

OPEN(UNlT=LLM,STATUS='OLO' ) 

OPEN (UN I T=LLN, FORMS' UNFORMATTED • , ST ATUS= * NEW ' ) 

KROrt IS THE NEW NUMB. OF ROWS ,KJ IS THE OLO 

NCOUN COUNTS NODAL PTS IN EPSA, LLID COUNTS FOR ACTUAL MOOEL 
NUMR COUNTS NOOAL PTS IN EPSA 

FIRST LOOP TO GET MAX. OEFORMATION AND NOOE NUMBER 

00 200 J=l,KROW 



SUBROUTINE NEUOISP(DEFL, VELO) 

OIMENSION A ( 1 ) 

COMMON I A ( 1 ) 

EQUIVALENCE ( I A ( 1 ) , A ( 1 > ) 

DIMENSION DEFLC 1 ) »VELOm 

COMMON /CPARA/ NSTEPS , NBEG , NENO , NSHEET , N2B0 , N3BD1 , 

1 N 3802 , INTRVL, OELT , NHTOT, HJOIN, NRELAX, ALPHA 
S LEN S8X 

COMMON / SSIZE / IBGOI , NO J , NEL TO T , N 1 80 , NLOA 0 , NBREC T , 

1 NBGUAD. ISHEET, NPRPTS, NSPTS, NSTRPTS, NVPTS, NHPTS, 

2 NSSTrP, NNJ, NNTOT , LGOSP, LIQUO, L8C ALC 

COMMON /STAB / IBT(1),JSSIZE,JSPAR,JVEL0,JSTRE,JXMAS,JIELM,JBMAT, 

1 JL1B0, JLOOP, JPRET, JLHI3, JSTRN, JFORC , JXLOC, JNQI , JNNI , JNQ8EG, 

$ JLSIOE/ JIELMCL.JSTIF, JDEFL.JFORLG, 

2 JIFPAR, JFLPAR, JXCOORO, JYCOORO.JOELTAX, JOEL T AY, JVM A, JSEFX, 

3 JFLUFR, JPRINC, JVELRAD, JGENFR , JPRES , JCSEP 

COMMON /CIO/ NIN,NOUT,NTHIST,NCORT,MCORT,NTPLOT,NVMA, ITITLE(20) 

COMMON/FRANK/NFNTOT,LLU,LLN,LLS,NOPLTS,NSPLTS 

COMMON/ T I TRE/N TITLE (80) 

COMMON/COELT/ISTEP, TIME 
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NCOjN=NCDUN+ I A ( JNNI + J-l ) 

LlPT=lA(JNNI+J)-2 

C 

C LLP T IS NUMB. OF POINTS (FOR P A T R A N ) IN EACH ROW 

C *E TAKE THE RIGID BODY MOTION OUT 8Y SUBTRACTING THE V AND W 

c displacements AT ENDBLOCK AT EACH NODAL POINT 
C NC OUN 1 2 IS THE POINT NB. (EPSA) FOR END9LDCK 

C 

R8Y=0EFL(3*(NCDUN+2)-l ) 

RBZ=DEFL(3*(NCOUN+2) ) 

IF (LIQUO.EQ.O) THEN 
RBY = 0 . 

RBZ=0. 

END I F 

DO 100 L= l r LLPT 
LL I D=LL ID* 1 
NUMR=NCOUN+L+ 1 

DX=DEFL(3*NUMR-2) 

DY=DEFL(3*NUMR-1 )-R8Y 
OZ=OEFL(3**NUMR)-RBZ 
C 

C IS SHELL CURVED, CHANGE COORDINATES 

C 

IF( (A ( JVELD-1 ) .NE.O. ) .DR. (A( JVELO-2) .NE.O.) ) THEN 
CALL CHCDDRD(NUMR,DX,DY,DZ) 

ENDIF 

00 = AM AX 1 (ABS(DX) , ABS(DY) ,ABS(DZ)) 

IF (ABS (DEFMAX) .LT.OO) THEN 

NOMA X =LL I 0 

OEFmax=OD 

IFCABSC AMIN1 (DX,DY,DZ) ) .EQ.DD) THEN 
OEFMAX=-OEFMAX 
END I F 
C 

C WE CHECKED IF DEFMAX WAS NEG. 

C 

END IF 

100 CONTINUE 

200 CONTINUE 

C OK FOR FIRST CARO 

WRITE(LLN) (N TITLE (I) , 1=1 ,80) , NFNTDT , MA XNDO , DEFMAX , NDMAX , 

1 NWIDTH 

wpi teCLLM, 90) T I TLE , NFN TOT , MA XNOO , DEFM AX , NOMAX , NW IOTH 

90 FORMATC A,2I8,E16.9,2I8) 

95 FORMA T (20A 1 ,218, E 16.9,21 8) 

WRITE(LLN) (NSU3T1 (I) , 1=1,80) 

WRITE(LLN) CNSUBT2(I), 1=1,80) 

WRITE(LLM,91 ) 

91 FORMATC’PINE* ) 

99 FORMAT(80A1) 

LL I D = 0 
N CDUN=0 
C 

C SECOND LOOP TO PICK UP DEFLECTION AT EACH NODE (OF PATRAN MODEL) 

C 

OD 400 J=1 , KROW 

NCD'JN = NCOUN + IA( JNNI f J-l ) 

LLPT=IA( JNNI+J)-2 
C 

C TAKE OUT RIGID BODY MOTION, R6X , RBZ DISPLACEMENTS AT END BLOCKS 

C AStJMEO TO REPRESENT RIGID BODY m 0T i 0 n S 
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Q8Y=0EFL(3* (NC0U.Nt2)-l ) 

RBZ=0EFL(3*(NC0lJN*2)) 

IF (LI QUO .EQ . 0 ) THEN 
RBY = 0 . 

RBZ=0. 

END I F 

C IF NO FLU 1 0 OPTION 00 NOT SUBTACT RIGIO BODY CONTRIBUTION 
00 300 L=1,LLPT 
LL 1 0 = LL 1 0 + 1 
NUMR=NC0UN+L+1 
0X=0EFL(3*NUMR-2) 

DY=DEFL( 3*NUMR-1 )-R9Y 
OZ=OEFL(3*NUMR)-R8Z 
OZLOC =0 Z 

VELZ=VEL0(3*NUMR) 

IF( ( A( JVELO-1 ) .NE.O. > .OR. ( A< JVELO-2) .NE.O. ) ) THEN 
CALL CHCOORO(NUMR,OX,OY,OZ) 

ENOIF 

rtRITE(LLN) LLIO,OX,DY,DZ,OZLOC, VELZ 
*IRITE(LLM,92) LLIO,OX,DY,DZ,OZLOC,VELZ 
300 CONTINUE 

aoo CONTINUE 

C CLOSE FILE OPENEO ON UNIT LLN 
CALL CLOSE ( LLN ) 

92 FORMAT (I8,5E16.9) 

RETURN 

ENO 

C 

C 

c 



SUBROUTINE CHCOORO ( N , X , Y , Z ) REPORT 

0 1 MENS ION A ( 1 ) * BLANC 

COMMON I A ( 1 ) BLANC 

EQUIVALENCE (IA(1), A ( 1 ) ) BLANC 

COMMON /CPARA/ NS T EPS , N8EG , NENO , NSHEET , N290 , N3B0 1 , CPARA 

1 N3802, INTRVL, DEL T , NHTOT, NJOIN, NRELAX, ALPHA CPARA 

$ LEN SB X MSP AR A 

COMMON / SSIZE / IBGCl ) , NO J , NEL TO T , N 1 BO , NLO AO , NBREC T , SSIZE 

1 NBQUAD, ISHEET, NPRPTS, NSPTS, NSTRPTS, NVPTS, NHPTS, SSIZE 

2 NSSTYP, MNJ , NNTOT, LGOSP, LIO'JD, LBC ALC SSIZE 

C REPORT 

COMMON /STA8 / I BT ( 1 ) , JSS IZE , JSP A R , J V ELO , JS TRE , JX M A S , J I ELM , JBM A T , STAB 

1 JL1B0, JLOOP, JPRET, JLHIS, JSTRN, JFORC, JXLOC, JNQI , JNNI , JNQBEG, STAB 

$ JLSIDE, JIELMCL, JSTIF, JDEFL, JFORLG, STAB 

2 JIFPAR, JFLPAR, JXCOORO, JY COORD, JOELTAX, JOELTAY, JVM a , JSEFX, STAB 

3 JFLUFR , JPRINC, JVELRA 0 , JGENFR , JPRES, JCSEP STAB 



COMMON/ CIO/ N I N, NOU T, N TH I 3 T, NIC OR T, MCO RT, NT PLOT, NVM A, I TITLE (20) 

common/frank/nfntot,llu,lln,lls,noplts,nsplts 

0 I MENS I ON XMAT(3,3) 

00 10 1=1,3 
00 10 J= t , 3 
XMAT(I, J)=0. 

10 CONTINUE 
C 

C THIS ROUTINE CHANGES THE OISPL. OF NOOAL PT. N ( FOR EPSA) 

C IN A GLOBAL RECTANGULAR SYSTEM 

C 

C IF CURVATURE IN Y OIR. IS NON ZERO /CALCULATE 

C THE ROTATION MATRIX AT EACH POINT 

C 



2 

2 

3 

a 
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3 
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IF (A(JvELO-l J.NE.O.) THEN 
PC3URY=1 . / A ( JVELO- 1 ) 

THET A1=A ( JXL0C*2*N-1) /RCOURY 
XM A T ( 1 , 1 ) = 1 . 

XMAT(2,2)=C0S(THETA1 ) 

XMAT ( 2 / 3)=SIN(THETA1 ) 

XMAT (3/3)-XMAT (2/2) 

XMAT(3r2)=-l . *XMA T ( 2 / 3 ) 

CALL PROD(X M AT / X / Y / Z ) 

END I F 
C 

C IF CURVATURE IN X DIRECTION IS NON ZERO: 

C 

IF ( A( JVELO-2) .NE.O. ) THEN 
DO 20 1=1/3 
DO 20 J= 1 / 3 

XMAT (I / J ) =0 . 

20 CONTINUE 

RC0UPX=1 ./A( JVELO-2) 

THETA2=A(JXL0C+2*N-2)/RC0URX 
XMAT (2/2) = l . 

XMAT (1 / 1 )=C0S(THETA2) 

XMAT(3/3)=XMAT(1/ 1) 

XMAT ( 3 / 1 )=-l .*SIN(THETA2) 

XMAT C 1 / 3)=SIN(THET A2) 

CALL PROD(XMAT,X/ Y,Z) 

END IF 
RETURN 
END 
C 
C 

SUBROUTINE PROD (X^AT/X/Y/Z) 

DIMENSION XMAT (3/ 3) 

C 

C THIS ROUTINE DOES THE MATRIX PRODUCT XMAT*(X/Y/Z) 

c 

X1=X 

Y 1 = Y 

Z 1 =z 

XsXMAT ( 1 , 1 ) *X1 tXMAT (1 ,2) *Y 1 tXMAT ( 1 , 3) *Z1 

Y = XMAT ( 2 / l)*Xl *XMAT(2,2)*Y1+XMAT(2,3) *Z1 
Z=XMAT(3, 1 )*Xl*XMAT(3,2)*YlfXMAT(3, 3)*Z1 
RETURN 

END 

C 

C 



CPARA 2 
CPARA 3 
MSP AR A 3 



C CHECKS IF OUTPUT FOR PATRAN IS REQUESTED BY USER IF YES CALL 

C NEUDISP 

C 

IF ( T IME • EQ . DEL T ) KDISP=1 
TIMEI=FLOAT(NSTEPS) /FLOAT (NDPLTS) 

TI^EC=KDISP*TIMEI 
ITIME=JNINT CTIMEC) 



SUBROUTINE PLOTD ISP ( DEFL / VELO ) 

DIMENSION DEFL ( 1 ) / VELO C 1 ) 

COMMON/CDELT/ISTEP/ TIME 
COMMON /CPARA/ NSTEPS , N9EG / MEND , 

1 N3BD2 * INTRVL, DELT/ NHTOT , NJOIN, 

$ LEN SB X 

COMMON/FRANK/NFNTOT ,LLU/LLN,LLS/ NDPLTS /NSPLTS 



NSHEET / N2BD 
NRELAX, ALPHA 



N3BD1 
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TIMEC=TIME/(ITIME*DELT) 

ENOTIMsnSTEPS*DELT 

TEST=ABS(TIMEC-JNINT(TIMECn 

IF( (TEST. LE. 0.0 0001 ) ,0R. (TIME.EQ.ENOTIM ) ) THEN 
CALL NEUDISP(0EFL, VELO) 

K0ISP=KDISP+1 

end i f 

RETURN 

ENO 

C 

SUBROUTINE PLOT STRE ( STRE ) 

DIMENSION STRE ( n 
COMMON/COELT/ISTEPr TIME 

COMMON /CPARA/ NSTEPS , NBEG , NEND , NSHEET , N2B0 , N3BD1 , CPAR A 

1 N3BD2, INTRVL, OELT , NHTOT , NJOIN, NREL AX , ALPHA CPARA 

$ LEN SB X MSP AR A 

common/frank/nfntot,llu,lln,lls,noplts,nsplts 

C CHECKS IF OUTPUT FOR PATRAN IS REQUESTEO BY USER IF YES CALL 

C NEUDISP 

C 

IF (TIME. EO. OELT) KSTRE=1 
T I ME I =FLO AT (NSTEPS) /FLO AT (NSPLTS) 

TIMEC=KSTRE*TIMEI 

ITIME=JNINT(TIMEC) 

TIMEC=TIME/( ITIME*OELT) 

ENDTIM=NSTEPS*DELT 
TEST=ABS(TIMEC-JNINT(TIMEC) ) 

IF C (TEST. LE. 0.0 0001 ) .OR. ( T I ME . EQ . ENOT I M ) ) THEN 
CALL NEUSTRE ( STRE ) 

KSTRE=KSTRE+ 1 
ENOIF 
RETURN 
END 
C 
C 



SUBROUTINE NEUSTRE ( STRE ) REPORT 

OIMENSION A ( 1 ) BLANC 

COMMON I A ( 1 *) BLANC 

EQUIVALENCE (IA(l)r A ( 1 ) ) BLANC 

OIMENSION OISP( 10) , STRA ( 1 0 ) , DSTRE ( 1 0 ) , OS ( 5 ) , VEL ( 2 ) 

OIMENSION STPE(l) 

COMMON /CPARA/ NSTEPS , NBEG , NEND , NSHEET , N2B0 , N3B01 , CPARA 

1 N3BD2, INTRVL, DELT, NHTOT, NJOIN, NREL AX , ALPHA CPARA 

$ LEN SBX MSP AR A 

COMMON / SSI2E / IBG( 1 ) ,NQJ,NELT0T,N1B0,NL0A0,NBRECT, SSIZE 

1 NBQU AO , ISHEET, NPRPTS, NSPTS, NSTRPTS, NVPTS, NHPTS, SSIZE 

2 NSSTYP, NNJ, NN TOT , LGOSP, LIQUO, LBCALC SSIZE 

C REPORT 

COMMON /STAB / I B T ( 1) , J SS I ZE , JSP AR , J V EL 0 , JS TRE , J XM A S , J I ELM , JBM AT , STAB 

1 JL1BD, JLOOP, JPRET, JLHIS, JSTRN, JFORC, JXLOC, JNQI , JNNI , JNQBEG, STAB 

$ JLSIOE, JIELMCL, JSTIF, JOEFL, JFORLG, STAB 

2 JIFPAR, JFLPAR, JXCOORO, JYCOORU, JOELTAX, JOELTAY, JVM A, JSEFX, STAB 

3 JFLUFR, JPRINCr JVELRAO , JGENFR , JPRES , JCSEP STAB 

COMMON/CIO/ NIN,NOUT,NTHIST,NCORT,MCORT,NTPLOT,NVMA, I T I TLE (20 ) 

COMMON/ FRANK/NFNTOT,LLU,LLN,LLS,NDPLTS, NSPLTS 

COMMON/TI TRE/NT ITLE (80 ) 

COMMON/COELT/ISTEP, TIME 



2 

3 

3 

2 
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COMMON/ SPAR A /HEG( l ) ,E,GNU,RHQ,SGYLO, THICK ,CURV (2) 
INTEGER NSU9T1 (80) ,NSU9T2(80) 

NWlDTH=31 
LLS= l 6 
LLZ=l5 

00 10 J=l,80 
10 NSUB T 2 ( J ) -0 

C TAKE CARE OF The 2N0 TITLE , IN0ICAT. TIME STEP 
C 

CL0S£(UNIT=LLZ) 

OPEN(UNIT=LLZ/STATUS=' NEW* ) 

WRITE (LLZ t * ) TIME 
CLOSE(UNIT=LLZ) 

RE AO ( LLZ #801) (NSU8TI (I) ,1 = 1,80) 

CLOSE(UNlT=LLZ) 

OPEN(UNIT=LLZ#STATUS=*OLO* ) 

801 FORMAT (80A1 ) 

OPEN(UNIT=LLS, FORMs 'UNFORMATTED ' , S T A TUS = ' NEW 1 ) 



WE ARE GOING TO DISPLAY THE FIRST THREE RECOROS (TITLES) 

WRITE (LLS) (NTITLE(I),t=l,80),NWIOTH 
WRITE (LLS) (NSUBTl (I) , 1 = 1,80) 

WRITE (LLS) (NSUBT2(I) , 1 = 1,80) 

WE ARE' GOING TO SCAN ALL ELEMENTS AS IN SHCONN, 

GET THE STRESSES PERFORM ROTATION IF NECESSARY 
LLEL: NUM. OF ROWS ELTS; LLIO : ELT NUMBER 
NPREV: NOOE NUM . ( FDR EPSA) OF LAST POINT OF ROW JUST 
BELOW THE ELEMENT ROW K; LLROW : NUM. OF ELTS IN EACH ROW 

NPREV=0 

IA ELTS + l (FOR NOOES) +2 ARTIFICIAL NOOES FOR NPREV 

NSHAPE=4 
LLEL=JNNI-JNQI 
00 200 Ks l , LLEL 

LLR0W=IA(JNQI+K-1) 

NPREV=NPREV*LLROW+3 

00 100 J=t, LLROW 
LLI0=IA( JNOBEG+K-l) +J-1 
OS ( l ) =STRE (9* (LLIO-1 ) +1 ) 

DS(2)=STRE(9*(LLIO-l)+2) 

OS ( 3) =0 . 

0S(4,5) WILL BE THE STRESSES IN LOCAL COORD. 

0S(4)=0S(l) 

OS ( 5 ) =OS ( 2 ) 

VON MISES AT TOP ANO BOTTOM OF SHELL: 

T0XY=STRE(9*(LLID-1 )+3) 

0SX=STRE(9*(LLI0-l)+a) 

DSY=STRE(9*(LLlD-l)+5) 

0SX=0SX*6. /THICK 
0SY=0SY*6. /THICK 
DSX l =OS ( l )+OSX 
OS Y l =OS ( 2 ) +0S Y 
DSX2=0S(l)-0SX 
OSY2=OS(2)-OSY 

VON MISES CALCULATIONS MODIFIED AFTERWARDS FOR ME MB. +8ENOING STRESS 
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NO SHEAR AT TOP OF SHELL , GET STRESSES AT TOP ANO BOTTOM AMO TAKE 
MAX. VALUE 

TOP: 

SIGMlrosxl 
S I GM2-OS Y 1 

V0NM1=SQRTCSIGM1 *SIGM1-SIGM1 *SIGM2+5IGM2*SIGM2) 
BOTTOM: 

SIGMl=OSX2 

SIGM2=0SY2 

VONM2=SQRT (SIGM1 *SIGM1-SIGM1 * S I GM2 + S I GM2 * S I GM2 ) 

VONM~ AM A X 1 (V0NM1 , VONM2) 

O.K. FOR VON MISES (VONM) 

IS THE SHELL CURVEO, NLELI 1 S ARE NOOES SURROUNOING ELT. (EPSA) 

IF ( (ACJVELO-l ) .NE.O.) .OR . ( A ( JVELO-2 ) .NE.O. ) ) THEN 

NLELl=J+NPREV+i 

NLEL2=J>NPREV+2 

NLEL3=J + H-NPREV+LLROW+H-2 

NLEL<i=J+NPREV+LLROW+ 1+2 

X=(A(JXLOC+2*NLELl-2)+A(JXLOC+2*NLEL2-2) )/u. 

X=X+(A (JXLOC +2*NLEL3-2)+A(JXLOC+2*NLELU-2) ) /U. 

Y=(A( JXLOC *2*NLEL1-1 ) + A ( JXLOC +2+NLEL2- 1 ) ) /4. 

Y=Y+( A (JXLOC t2*NLEL 3-1 )*A( JXLOC +2*NLELU-1 ) )/4. 

CALL CHELEMCX, Y,0S) 

ENOIF 

REAOY TO DISPLAY RECORO 



C 

100 

200 

80 



WRITE(LLS) LLIO ,NSHAP£, (OISP ( I ) , 1=1 , 1 0) , 

1 (STRA(I) , 1st, 10) ,OSTRE(l) ,OS(5) ,0S(6) , 

2 OS T RE ( 4 ) , C OS ( I ) , 1 = 1/3)# 

3 (OSTRE(I) ,1=8,10) , VONM 

WRITE (15,30) LLIO , NSHAPE , (OS ( I ) , 1*1 , 5) ,V0NM 
CONTINUE 
CONTINUE 
CALL CLOSE (LLS ) 

FORMAT(2ia,6El*3.8) 

RETURN 

ENO 



C 



SUBROUTINE CHELEMCX, Y,0S) 

OIMENSION A ( 1 ) 

COMMON I A ( 1 ) 

EQUIVALENCE ( I A ( 1 ) , A c X ) ) 

COMMON /CPARA/ NSTEPS , N8EG , NENO , NSHEET , N20O , N3B01 , 
l N3B02, INTRVL, OELT , NHTOT , NJOIN, NRELAX, ALPHA 
$ LEN S3 X 

COMMON / SSIZE / I9G( 1 ) ,NQJ ,NELTOT ,N190,NL0A0,N6RECT , 

1 NBQUAO, ISHEET, NPRPTS, NSPTS, NSTRPTS , NVPTS, NHPTS, 

2 NSSTYP, NNJ, NNTOT, LGOSP, LIQUO, LBCALC 

COMMON /STAB / IBT(l),JSSIZE,JSPAR,JVELO,JSTRE,JXMAS,JIELM,JBMAT, 

1 JL1B0, JLOOP, JPRET , JLHIS, JSTRN, JFORC, JXLOC, JNOI, JNNI, JNQBEG, 

S JLSIOE, JIELMCL, JSTIF, JOEFL, JFORLG, 

2 JIFPAR, JFLPAR, JXCOORO, JYCOORO, JOELTAX, JOELTAY, JVM A, JSEFX, 

3 JFLUFR, JPRINC, JVELRAO, JGENFR , JPRES , JCSEP 
OIMENSION XMAT(3,3) 

OIMENSION OS ( 3 ) 



REPORT 2 
BLANC 2 
BLANC 3 
BLANC a 
CPARA 2 
CPARA 3 
MSP AR A 3 
SSIZE 2 
SSIZE 3 
SSIZE a 
REPORT 8 
STAB 2 
STAB 3 
STA8 a 
STAB 5 
STAB 6 
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oo to 1=1,3 
00 10 j=i , 3 
XMAT ( i , J) =o. 

10 CONTINUE 

THIS ROUTINE CHANGES THE STRESSE IN AN ELEMENT INTO A 
IN A GLOBAL RECTANGULAR SYSTEM .GIVEN LOCAL CQORO. OF 
CENTRO 1 0 X , Y 



IF CURVATURE IN Y OIR. IS NON ZERO 
THE ROTATION MATRIX AT EACH POINT 
IF (A(JVELO-l).NE.O.) THEN 
RCOtJRY = 1 ./A(JVELO-t) 

THETA 1=Y/RC0URY 
XMAT ( l , 1 )=1 . 

XMAT(2,2)=C0S(Th£TAI ) 

XMAT (2,3)=SIN(TH£TA1 ) 



XMAT (3,3)=XMAT (2,2) 
XMAT(3,2)=-1.*XMAT(2,3) 

CALL PROO(XMAT,OS( 1 ) ,0SC2) , OS C 3 } ) 
ENOIF 
C 

C IF CURVATURE IN X OIRECTION IS NON ZERO: 
IF ( A( JVELO-2) .NE.O. ) THEN 
00 20 1=1,3 
00 20 J= 1 , 3 
XMAT ( I , J ) =0 . 

20 CONTINUE 

RCOURXs 1 ./A( JVELO-2) 
THETA2=X/RC0URX 
XMAT (2,2)=1 . 

XMAT ( 1 , t )=C0S(THETA2) 

XMAT (3, 3)=XMAT ( 1,1) 

XMAT ( 3 , 1)=-1.*SIN(THETA2) 

XMAT (1,3)=SIN (THETAS) 

CALL PROD (XMAT ,0S( 1 ) ,0S(2) , 0S(3) ) 
ENOIF 
RETURN 
ENO 



.calculate 
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